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The  impact  of  parameterized  topographic  internal  lee  wave  drag  on  the  input  and  output  terms  in  the 
total  mechanical  energy  budget  of  a  hybrid  coordinate  high-resolution  global  ocean  general  circulation 
model  forced  by  winds  and  air-sea  buoyancy  fluxes  is  examined  here.  Wave  drag,  which  parameterizes 
the  generation  of  internal  lee  waves  arising  from  geostrophic  flow  impinging  upon  rough  topography, 
is  included  in  the  prognostic  model,  ensuring  that  abyssal  currents  and  stratification  in  the  model  are 
affected  by  the  wave  drag. 

An  inline  mechanical  (kinetic  plus  gravitational  potential)  energy  budget  including  four  dissipative 
terms  (parameterized  topographic  internal  lee  wave  drag,  quadratic  bottom  boundary  layer  drag,  vertical 
eddy  viscosity,  and  horizontal  eddy  viscosity)  demonstrates  that  wave  drag  dissipates  less  energy  in  the 
model  than  a  diagnostic  (offline)  estimate  would  suggest,  due  to  reductions  in  both  the  abyssal  currents 
and  stratification.  The  equator  experiences  the  largest  reduction  in  energy  dissipation  associated  with 
wave  drag  in  inline  versus  offline  estimates.  Quadratic  bottom  drag  is  the  energy  sink  most  affected  glob¬ 
ally  by  the  presence  of  wave  drag  in  the  model;  other  energy  sinks  are  substantially  affected  locally,  but 
not  in  their  global  integrals.  It  is  suggested  that  wave  drag  cannot  be  mimicked  by  artificially  increasing 
the  quadratic  bottom  drag  because  the  energy  dissipation  rates  associated  with  bottom  drag  are  not  spa¬ 
tially  correlated  with  those  associated  with  wave  drag  where  the  latter  are  small.  Additionally,  in  contrast 
to  bottom  drag,  wave  drag  is  a  non-local  energy  sink. 

All  four  aforementioned  dissipative  terms  contribute  substantially  to  the  total  energy  dissipation  rate 
of  about  one  terawatt.  The  partial  time  derivative  of  potential  energy  (non-zero  since  the  isopycnal 
depths  have  a  long  adjustment  time),  the  surface  advective  fluxes  of  potential  energy,  the  rate  of  change 
of  potential  energy  due  to  diffusive  mass  fluxes,  and  the  conversion  between  internal  energy  and  poten¬ 
tial  energy  also  play  a  non-negligible  role  in  the  total  mechanical  energy  budget.  Reasons  for  the  <10% 
total  mechanical  energy  budget  imbalance  are  discussed. 

©  2013  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

In  this  paper,  we  investigate  the  impact  of  parameterized  topo¬ 
graphic  internal  lee  wave  drag  on  the  input  and  output  terms  in  the 
total  mechanical  energy  budget  of  a  high-resolution  (“eddying”) 
global  ocean  model.  We  are  motivated  by  the  potentially  important 
role  of  topographic  internal  lee  wave  drag  in  mixing  the  abyssal 
ocean.  In  recent  years,  there  has  been  great  interest  in  the  ocean 
energy  budget,  largely  because  the  mixing  associated  with  energy 
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dissipation  is  thought  to  exert  an  important  control  on  the  large- 
scale  circulation.  Munk  and  Wunsch  (1998)  and  St.  Laurent  and 
Simmons  (2006)  have  suggested  that  about  2-3  TW  of  mixing  en¬ 
ergy  is  required  to  raise  diffusivity  enough  to  maintain  the  abyssal 
stratification  in  the  presence  of  the  30  Sverdrups  ( Sv  =  10s  m3  s_1) 
of  deep  water  formation.  However,  Webb  and  Suginohara  (2001) 
have  suggested  that  maintaining  9  Sv  of  Ekman  suction  in  the 
Southern  Ocean  while  vertically  mixing  17Sv  of  North  Atlantic 
Deep  Water  would  reduce  the  required  energy  dissipation  rate  in 
the  abyssal  ocean  to  as  little  as  0.6  TW. 

Recently,  intense  research  interest  has  focused  on  the  sources 
and  sinks  of  mixing  energy  in  the  ocean.  Almost  all  of  the  60-68 
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terrawatts  (TW)  of  wind  power  put  into  the  surface  waves  is  dissi¬ 
pated  near  the  surface  (Wang  and  Huang,  2004;  Ferrari  and 
Wunsch,  2010)  and  most  of  this  wind  power  input  is  thought  to  en¬ 
hance  vertical  shear  of  the  mean  currents  (Large  et  al.,  1994)  and 
vertical  mixing  (Wang  et  al.,  2010;  Shu  et  al.,  2011;  Qiao  and 
Huang,  2012)  in  the  upper  ocean.  About  0.4  TW  of  wind  power  is 
put  into  near-inertial  motions  in  the  mixed  layer  (Watanabe  and 
Hibiya,  1239;  Alford,  2003;  Furuichi  et  al.,  2008).  Early  estimates 
(Wunsch,  1998;  Scott  and  Xu,  2009)  of  the  wind  power  put  into 
the  ocean  general  circulation,  which  includes  low-frequency  cur¬ 
rents  and  mesoscale  eddies,  were  found  to  be  about  0.9  TW.  How¬ 
ever,  Zhai  et  al.  (2012)  suggest  that  only  about  0.5  TW  of  wind 
power  is  put  into  the  general  circulation  because  higher  frequency 
wind  variability  generates  shear,  mixing,  and  near-inertial  waves 
in  the  surface  layer  rather  than  deep  ocean  mixing. 

The  energy  sinks  of  the  low-frequency  eddying  oceanic  general 
circulation  are  just  beginning  to  be  estimated  on  a  global  scale. 
Topographic  internal  lee  wave  drag  is  a  potentially  critical  compo¬ 
nent  of  the  mechanical  energy  budget.  Naveira-Garabato  et  al. 
(2004),  Marshall  and  Naveira-Garabato  (2008),  and  Nikurashin 
(2008)  suggested  that  energy  is  transferred  to  internal  lee  waves 
when  geostrophic  flow  impinges  upon  rough  topographic  features 
and  is  eventually  dissipated,  especially  in  the  Southern  Ocean 
where  geostrophic  flows  are  strong  and  the  bottom  is  rough.  This 
energy  dissipation  mechanism,  which  is  the  main  focus  of  the  pres¬ 
ent  study,  will  be  referred  to  simply  as  "wave  drag”  hereafter. 

Other  postulated  energy  dissipation  mechanisms  for  the  eddy¬ 
ing  general  circulation  include  quadratic  bottom  boundary  layer 
drag  (hereafter,  “bottom  drag”;  Sen  et  al.,  2008;  Arbic  et  al.,  2009; 
and  references  therein);  energy  scattering  into  high-wavenumber 
vertical  modes  (Zhai  et  al.,  2010;  Saenko  et  al.,  2012);  and  cata¬ 
lyzed  energy  exchanges  via  inviscid  balanced  flow-boundary  inter¬ 
action  (Dewar  and  Hogg,  2010;  Dewar  et  al.,  2011).  In  ocean 
models,  energy  is  also  dissipated  by  the  vertical  eddy  viscosity 
(Large  et  al.,  1994)  and  horizontal  eddy  viscosity  (Smagorinsky, 
1993)  that  must  be  employed  to  make  up  for  the  lack  of  resolved 
small-scale  processes.  Vertical  eddy  viscosity  (“vertical  viscosity" 
hereafter)  represents  processes  associated  with  vertical  shear 
instabilities.  Horizontal  eddy  viscosity  (“horizontal  viscosity”  here¬ 
after)  is  meant  to  represent  processes  that  can  remove  vorticity 
and  momentum  at  the  boundaries  of  ocean  basins  (Fox-Kemper 
and  Pedlosky,  2004).  Arguably,  horizontal  viscosity  also  very 
roughly  represents  small-scale  processes  -  for  instance,  energy 
transfer  from  mesoscale  eddies  to  either  internal  waves  (Polzin, 
2008)  or  submesoscale  eddies  (Muller  et  al.,  2005)  -  which  are 
not  explicitly  resolved  by  any  existing  numerical  global  ocean 
model. 

Here,  we  quantify  the  relative  amounts  of  energy  dissipation  of 
low-frequency  flow  in  an  eddying  ocean  model  due  to  bottom  drag, 
wave  drag,  horizontal  viscosity,  and  vertical  viscosity.  Previous 
estimates  suggest  that  bottom  drag  and  wave  drag  both  contribute 
substantially  to  the  energy  budget  of  low-frequency  flows.  For 
example,  Sen  et  al.  (2008)1,  Wright  et  al.  (2012), 2  and  Arbic  et  al. 
(2009 )3  have  argued  that  bottom  drag  dissipates  at  least  0.2  TW  of 
low-frequency  mechanical  energy.  Arbic  and  Flierl  (2004)  and 
Wright  et  al.  (2013)  further  argued  that  some  of  the  energy  dissipa¬ 
tion  that  is  typically  attributed  to  bottom  drag  in  both  models  and 
observations  should  actually  be  attributed  to  wave  drag.  Nikurashin 
and  Ferrari  (201 1 )  estimated  the  rate  of  energy  dissipation  by  break¬ 
ing  lee  waves  to  be  about  0.2  TW  by  assuming  that  this  rate  is  a  frac¬ 


1  They  made  use  of  the  Deep  Water  Archive  and  Buoy  Group  Archive,  (http:// 
cmdac.oce.orst.edu/cds.html  or  http://cmrecords.net) 

2  They  made  use  of  the  more  extensive  Global  Multi-Archive  Current  Meter 
Database  (http://stockage.univ-brest.fr/scott/GMACMD/updates.html). 

3  They  utilized  multiple  eddying  ocean  models. 


tion  of  the  energy  conversion  rate  into  internal  lee  waves.  In 
contrast,  Scott  et  al.  (2011)  estimated  the  rate  of  energy  conversion 
into  internal  lee  waves  to  be  about  0.34-0.49  TW.  Both  the 
Nikurashin  and  Ferrari  (2011)  and  Scott  et  al.  (2011)  estimates  are 
based  on  bottom  stratification  fields  taken  from  observations  in 
concert  with  bottom  flows  in  global  ocean  models  which  did  not  uti¬ 
lize  wave  drag. 

A  key  contribution  of  the  present  manuscript  is  the  prognostic 
calculation  of  wave  drag  and  an  evaluation  of  the  alterations  in  ki¬ 
netic  energy,  stratification,  and  sources  and  sinks  of  the  total 
mechanical  energy  budget  due  to  wave  drag  implementation. 
Our  insertion  of  wave  drag  into  a  global  ocean  general  circulation 
model  is  motivated,  in  part,  by  demonstrations  (e.g.,  Jayne  and 
St.  Laurent,  2001),  that  the  energy  budget  and  accuracy  of  global 
forward  tide  models  are  impacted  to  first  order  by  wave  drag.  Jay¬ 
ne  and  St.  Laurent  (2001 )  found  that  roughness  sufficient  to  gener¬ 
ate  internal  lee  waves  occurs  quite  commonly  in  the  open  ocean. 
Because  bottom  drag  only  depends  on  the  bottom  velocities  and 
not  the  roughness,  a  model  simulation  that  includes  a  wave  drag 
parameterization  should  have  more  dissipation  in  open  ocean  re¬ 
gions  than  a  model  simulation  that  only  includes  a  bottom  drag 
parameterization.  In  addition,  there  is  ample  evidence  from  both 
observations  (Polzin  et  al.,  1997;  Naveira-Garabato  et  al.,  2004; 
St.  Laurent  et  al.,  2012)  and  very  high-resolution  numerical  ocean 
process  models  that  include  bottom  drag  and  wave  drag  (Nikura¬ 
shin  et  al.,  2013)  that  turbulent  mixing  is  enhanced  when  low-fre¬ 
quency  flows  encounter  rough  topography. 

While  global  ocean  general  circulation  models  tend  to  be  defi¬ 
cient  in  bottom  kinetic  energy  relative  to  current  mooring  observa¬ 
tions  (Scott  et  al.,  201 1 ),  the  correlation  between  vertical  profiles  of 
kinetic  energy  in  ocean  models  versus  observations  may  be  a  more 
important  statistic  to  improve.  This  is  because  each  ocean  model 
grid  point  represents  an  average  over  a  large  area,  thus  tending 
to  smooth  the  kinetic  energy  at  each  model  grid  point  relative  to 
the  points  at  which  current  meter  measurements  are  taken.  Ocean 
models’  simulated  kinetic  energy  increases  with  finer  model  reso¬ 
lutions  (Thoppil  et  al.,  2011).  Therefore,  in  ocean  model  simula¬ 
tions  without  wave  drag,  bottom  kinetic  energy  may  be  closer  to 
that  of  current  mooring  observations  than  ocean  model  simula¬ 
tions  with  wave  drag  for  the  wrong  reasons  (i.e„  inadequate  reso¬ 
lution  in  combination  with  a  lack  of  abyssal  drag  such  as  wave 
drag).  It  will  be  left  to  a  future  manuscript  to  discuss  whether 
the  correlation  between  the  kinetic  energy  profiles  in  ocean 
models  versus  current  meter  observations  is  improved  with  the 
addition  of  wave  drag. 

In  this  paper,  we  analyze  the  global  total  mechanical  energy 
budget  of  the  total  (mean  plus  eddy)  flow,  using  global  nominally 
1/12°  simulations  of  the  HYbrid  Coordinate  Ocean  Model  (HYCOM; 
http://www.hycom.org;  Bleck,  2002;  Chassignet  et  al.,  2003; 
Halliwell,  2004)  with  and  without  wave  drag.  We  will  not  analyze 
the  following:  (1)  the  eddy  kinetic  energy  budget,  as  was  done  by 
Treguier  (1992)  using  a  1/3°  x  2/5°  ocean  model  of  the  North 
Atlantic;  (2)  the  generation  and  conversion  rates  between  gravita¬ 
tional  potential  energy  and  kinetic  energy,  as  was  done  by  Oort 
et  al.  (1994)  using  observations;  (3)  the  mean  kinetic  energy  and 
gravitational  potential  energy,  as  was  done  by  Aiki  et  al.  (2011) 
using  a  1/10°  global  ocean  model;  (4)  the  Lorenz  oceanic  energy 
cycle,  as  was  done  by  von  Storch  et  al.  (2012)  using  a  1/10°  global 
ocean  model;  or  (5)  the  kinetic  plus  available  potential  energy  bud¬ 
get,  as  was  done  by  Hogg  et  al.  (2013)  using  an  idealized  1/4°  ocean 
model  that  mimicked  the  Atlantic  Ocean.  The  conversion  of  kinetic 
energy  to  potential  energy  involves  work  done  by  several  pro¬ 
cesses,  some  of  which  include  horizontal  pressure  gradients,  verti¬ 
cal  velocities  that  result  from  the  convergence  or  divergence  of 
both  the  barotropic  and  baroclinic  components  of  the  horizontal 
velocities,  and  Reynolds  stresses  that  are  mediated  by  eddy  kinetic 
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energy  (Aiki  et  al.,  201 1 ).  The  generation  of  mean  kinetic  plus  po¬ 
tential  energy  is  approximately  balanced  by  the  energy  dissipators, 
and  the  baroclinic  energy  conversion  pathway  is  facilitated  by  the 
conversion  of  mean  kinetic  energy  to  mean  potential  energy  with 
the  surface  buoyancy  fluxes  and  winds  (von  Storch  et  al.,  2012). 
Here,  we  look  at  the  kinetic  plus  gravitational  potential  (i.e., 
mechanical)  energy  budget. 

This  manuscript  is  organized  as  follows.  First,  we  describe  the 
two  HYCOM  simulations  (with  versus  without  wave  drag).  Then 
we  provide  a  brief  description  of  the  quadratic  bottom  boundary 
layer  drag  parameterization.  Next,  we  introduce  the  abyssal  hill 
rough  topography  spectra  of  Goff  and  Jordan  (1988),  Goff  and  Arbic 
(2010),  and  Goff  (2010),  which  are  used  in  the  wave  drag  schemes. 
We  then  describe  two  different  wave  drag  schemes  utilized  in  this 
paper.  Our  analysis  with  the  Bell  (1975,  B75  hereafter)  wave  drag 
scheme  uses  the  bottom  densities,  stratification,  and  velocities  rel¬ 
ative  to  wavevectors  in  the  topographic  field’s  spectral  domain  to 
infer  an  energy  conversion  rate  into  lee  waves  (assumed  here  to 
be  the  energy  dissipation  rate  associated  with  wave  drag).  Our 
analysis  with  the  Gamer  (2005;  G05  hereafter)  wave  drag  scheme 
uses  the  stratification,  velocities,  and  features  of  the  underlying 
topography.  Because  the  G05  scheme  can  be  implemented  in  the 
physical  domain,  it  is  faster  and  easier  to  use  inline  with  a  global 
ocean  model.  After  describing  the  wave  drag  schemes  and  their  re¬ 
quired  input  parameters,  we  perform  offline  estimates  of  energy 
dissipation  rates  associated  with  wave  drag,  as  in  Nikurashin  and 
Ferrari  (201 1 )  and  Scott  et  al.  (201 1 ).  The  offline  estimates  are  use¬ 
ful  because  they  are  much  less  time-consuming  to  compute  than 
inline  estimates.  The  offline  estimates  provide  a  feasible  testing 
ground  for  comparing  the  G05  and  B75  schemes,  and  for  assessing 
the  impact  of  using  bottom  stratification  fields  taken  from  HYCOM 
as  opposed  to  observations.  We  then  run  a  HYCOM  simulation  with 
the  G05  scheme  inserted,  to  find  inline  estimates  of  the  total 
mechanical  energy  budget  derived  here,  including  energy  dissipa¬ 
tion  rates  due  to  wave  drag,  bottom  drag,  vertical  viscosity,  and 
horizontal  viscosity.  We  will  compare  the  relative  magnitudes  of 
the  total  mechanical  energy  budget  terms,  and  we  will  evaluate 
how  well  they  balance. 


2.  Model  description 

The  1/12°  HYCOM  simulations  are  on  a  global  tripolar  Mercator 
grid  (Murray,  1996)  and  have  32  hybrid  layers  in  the  vertical  direc¬ 
tion.  In  the  open  ocean,  the  model  employs  z-level  coordinates  in 
the  mixed  layer  and  isopycnal  coordinates  below  the  mixed  layer. 
In  shallow  areas,  the  model  employs  terrain-following  coordinates. 
This  hybrid  choice  is  motivated  by  the  strengths  of  the  different 
systems  in  their  respective  regions  (Griffies  et  al.,  2000).  Virtual 
potential  density  is  used,  referenced  to  2000  m  depth.  Unlike  po¬ 
tential  density,  virtual  potential  density  allows  the  inclusion  of 
buoyancy  anomalies  caused  by  thermobaricity  by  accounting  for 
compressibility  when  accelerations  due  to  pressure  gradients  are 
computed.  As  described  in  Sun  et  al.  (1999),  the  advantages  of  a 
virtual  potential  density  include  a  simple  pressure-gradient  force 
expression  and  monotonicity  in  depth.  Virtual  density  has  the  dis¬ 
advantages  of  being  nonmaterial  in  adiabatic  flow,  and  of  following 
neutrally  buoyant  surfaces  no  better  than  potential  density  does.4 

We  now  briefly  discuss  the  implementation  of  eddy  viscosities 
in  HYCOM.  The  K-Profile  Parameterization  (KPP;  Large  et  al.,  1994) 
is  used  to  determine  the  vertical  viscosity  fields,  which  are  derived 
from  the  vertical  diffusivities  with  an  assumed  Prandtl  number.  For 


4  Hallberg  (2005)  demonstrated  that  the  choice  made  by  Sun  et  al.  (1999)  to  define 
virtual  potential  density  can  lead  to  accelerations  due  to  pressure  gradients  that  are 
numerically  unstable  in  regions  with  weak  stratification. 


background  mixing,  which  typically  is  used  in  deep  water,  a  Pra¬ 
ndtl  number  of  three  is  assumed  so  that  the  vertical  viscosity  is  a 
factor  of  three  larger  than  the  vertical  diffusivity.  For  shear  insta¬ 
bility  mixing,  which  typically  is  used  in  the  mixed  layer,  a  Prandtl 
number  of  one  is  assumed.  KPP  yields  relatively  strong  vertical 
mixing  in  the  mixed  layer,  with  a  smooth  transition  to  weaker  ver¬ 
tical  mixing  below.  The  horizontal  viscosity  includes  the  maximum 
of  a  Laplacian  and  a  Smagorinsky  (1993)  parameterization  with  an 
additional  biharmonic  term.  Horizontal  viscosity,  employed  along 
layers  in  HYCOM,  smoothes  out  subgrid-scale  noise.  However, 
increasing  horizontal  viscosity  generally  comes  at  the  expense  of 
model  accuracy  (Wallcraft  et  al.,  2005;  Jochum  et  al.,  2008).  Here, 
“horizontal”  means  along-isopycnal  below  the  mixed  layer  in  the 
open  ocean,  along  a  constant  depth  surface  within  the  mixed  layer 
in  the  open  ocean,  and  parallel  to  the  bathymetry  where  terrain¬ 
following  coordinates  are  employed.  Referring  to  the  horizontal 
viscosity  term  as  “lateral  viscosity”  may  be  more  appropriate, 
but  we  refer  to  this  term  as  horizontal  viscosity  to  be  consistent 
with  our  reference  to  vertical  viscosity,  which  may  be  more 
accurately  called  a  “diapycnal  viscosity”.  For  offline  analysis,  the 
viscosity  terms,  the  temperature  and  salinity  fields  used  to  com¬ 
pute  the  stratification,  and  the  kinetic  energy  all  must  be  vertically 
interpolated  from  hybrid  coordinate  space  to  z-levels.  Each  three- 
dimensional  variable  in  the  kth  layer  (k  =  1, . . .  ,32)  included  in  the 
output  file  from  our  simulations  has  been  weighted  by  the  thick¬ 
ness  of  the  kth  layer  at  each  baroclinic  time  step  and  was  divided 
by  the  average  thickness  over  the  present  day  of  the  simulation 
when  the  output  was  written  to  file.  The  cumulative  sums  of  the 
thicknesses  are  the  depths,  so  these  are  easily  interpolated  to  a  reg¬ 
ular  depth  grid,  which  is  the  grid  we  present  our  zonally  averaged 
results  on. 

The  global  model  simulations  were  spun-up  from  rest  using 
1.125°  x  1.125°  European  Centre  for  Medium-Range  Weather 
Forecasts  (ECMWF)  Re-Analysis  (ERA-40)  monthly  mean  thermal 
forcing  over  1978-2002  (Kallberg  et  al.,  2004).  In  order  to  supple¬ 
ment  the  climatological  wind  forcing  with  higher  frequencies, 
six-hourly  anomalies  with  respect  to  monthly  means  from  the 
2003  fields  of  the  Navy  Operational  Global  Atmospheric  Prediction 
System  (NOGAPS;  Rosmond  et  al.,  2002)  was  added  to  the  ERA-40 
climatological  wind  forcing.  The  six-hourly  output  of  2003  NO¬ 
GAPS  winds  were  cycled  through  every  model  year  in  this  way. 
In  total,  HYCOM  was  spun-up  for  thirteen  years  without  wave  drag. 
To  accelerate  the  model  spin-up  period,  the  assumed  background 
tidal  velocity  (see,  for  instance,  Willebrand  et  al.,  2001)  was  not 
constant  throughout  the  thirteen  years.  Starting  from  rest,  the 
background  tidal  velocity  was  5  cm  s~]  for  the  first  one  and  one- 
half  years,  2  cm  s_1  for  the  next  two  and  one-half  years,  and 
0  cm  s_1  for  the  last  nine  years.  Following  the  thirteen  year 
spin-up  phase  without  wave  drag,  an  additional  simulation  was 
performed  for  another  seven  years  with  an  inserted  wave  drag 
scheme  using  an  assumed  tidal  velocity  of  0  cm  s_1.  For  an  offline 
analysis  of  wave  drag  using  two  different  schemes,  we  used  five- 
day  averaged  outputs  from  one  additional  (fourteenth)  year  of 
the  simulation  without  wave  drag. 

Our  inline  total  mechanical  energy  budget  analysis  was  per¬ 
formed  on  this  HYCOM  simulation  using  the  final  year  of  the  sim¬ 
ulation  without  wave  drag  (i.e.,  on  the  fourteenth  year)  and  using 
the  final  year  of  the  simulation  with  wave  drag  (i.e.,  on  the  twen¬ 
tieth  year).  We  found  that  results  computed  from  the  nineteenth 
year  of  the  spin-up  were  nearly  equal  to  results  from  the  twentieth 
year.  This  suggests  that  transient  adjustments  due  to  the  introduc¬ 
tion  of  wave  drag,  or  to  the  changing  background  tidal  velocities 
during  spin-up,  are  generally  not  large  by  the  time  we  compute 
our  energy  budgets.  The  globally  integrated  mean  kinetic  energy, 
shown  in  Fig.  1,  demonstrates  that  the  simulation  without  wave 
drag  is  reasonably  well  spun-up  after  thirteen  years  and  that  the 
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Depth-integrated  area-averaged  kinetic  energy  [kg  s'2] 


Fig.  1.  The  global  mean  (depth-integrated,  area-averaged)  kinetic  energy  [kgs-2], 
f  dzf  dA{0.5(u2  +  is2)p}/  f  dA,  as  a  function  of  spin-up  year  of  the  simulation 
without  wave  drag  (years  5-13)  and  continuation  upon  the  addition  of  wave  drag 
(years  14-19).  Here,  p  is  the  density,  f  dz  indicates  an  integral  over  all  depths,  and 
f  dA  indicates  an  integral  over  all  horizontal  grid  locations. 

simulation  with  wave  drag  is  reasonably  well  spun-up  after  an 
additional  six  years.  An  additional  diagnostic,  the  area-weighted 
center  of  mass  of  the  ocean,  also  suggests  that  the  model  is  reason¬ 
ably  well  spun-up  (not  shown).  While  the  globally  integrated  ki¬ 
netic  energy  decreases  upon  insertion  of  wave  drag,  the  globally 
integrated  potential  energy  increases  upon  insertion  of  wave  drag. 
We  will  evaluate  the  time  derivatives  of  the  globally  integrated  ki¬ 
netic  energy  and  gravitational  potential  energy  to  address  how 
important  transient  effects  are  to  our  total  mechanical  energy 
budget. 

Observational  evidence  that  wave  drag  should  be  applied  in  re¬ 
gions  where  there  are  abyssal  hills  comes  from  St.  Laurent  et  al. 
(2012)  and  Waterman  et  al.  (2013).  A  critical  decision  regarding 
simulations  with  wave  drag  is  what  to  do  in  non-abyssal  hill  re¬ 
gions,  defined  here  as  the  locations  where  the  abyssal  hill  rough 
topography  spectra  of  Goff  and  Jordan  (1988)  are  not  defined.  Be¬ 
cause  the  abyssal  hill  rough  topography  spectra  (see  Section  3.2) 
utilized  here  are  strictly  valid  only  for  abyssal  hill  regions,  we  at¬ 
tempted  a  model  simulation  in  which  wave  drag  was  applied  only 
in  abyssal  hill  regions.  However,  we  found  that  this  simulation  was 
not  numerically  stable  because  adjacent  grid  points  attempt  to  ap¬ 
ply  drags  that  are  orders  of  magnitude  different  from  each  other. 
Complementary  evidence  that  the  addition  of  a  wave  drag  field 
with  sharp  spatial  gradients  can  lead  to  numerical  problems  was 
found  in  Arbic  et  al.  (2010),  who  applied  wave  drag  to  tides  embed¬ 
ded  in  an  eddying  HYCOM  simulation.  Note  the  spuriously  large 
velocities  in  the  Gulf  of  Mexico  and  other  regions  in  Fig.  2b  of  Arbic 
et  al.  (2010)  relative  to  their  Fig.  2a.  In  subsequent  versions  of  HY¬ 
COM  simulations  with  embedded  tides  (e.g.,  Shriver  et  al.,  2012),  it 
was  found  that  smoothing  the  wave  drag  fields  eliminated  such 
artifacts.  All  of  this  suggests  that  applying  wave  drag  in  non-abys- 
sal  hill  regions  leads  to  smoother  numerics  in  the  model. 

The  question  then  becomes:  is  it  physically  justifiable  to  apply 
wave  drag  in  regions  without  abyssal  hills?  In  other  words,  is 
roughness  in  non-abyssal  hill  regions  substantial  enough  to  gener¬ 
ate  lee  waves?  Observational  justification  for  applying  wave  drag 
in  non-abyssal  hill  regions  comes  from  Kunze  et  al.  (2012)  and 
Beaird  et  al.  (2012).  The  former  study  found  that  the  breaking  of 
internal  waves  on  ridges  and  in  canyons,  which  can  be  found  in 
non-abyssal  hill  regions,  can  be  a  significant  sink  of  energy  (an  or¬ 
der  of  magnitude  more  energy  dissipated  locally  than  in  the  open 
ocean)  and  an  enhancer  of  diffusivity  (two  to  three  times  higher  lo¬ 


cally  than  in  the  open  ocean).  A  microstructure  survey  of  the  Ice- 
land-Faroe  Ridge,  performed  in  the  latter  study,  found  that  the 
turbulent  kinetic  energy  dissipation  was  enhanced  downstream 
of  the  Faroe  Bank  Channel  sills  and  in  the  overflow  region  along 
the  shelf  break  of  Iceland.  Neither  of  these  observational  studies 
found  that  lee  waves  are  necessarily  generated  in  their  respective 
regions,  but  they  do  suggest  that  energy  dissipation  is  enhanced 
in  non-abyssal  hill  regions  with  rough  topography.  Because  of 
the  numerical  problems  we  encountered  by  masking  out  wave 
drag  in  non-abyssal  hill  regions  and  the  observational  evidence 
that  there  is  enhanced  dissipation  in  non-abyssal  hill  regions,  here 
we  choose  to  extend  the  wave  drag  parameterization  into  non- 
abyssal  hill  regions. 

Strictly  speaking,  our  wave  drag  parameterizations  apply  only 
to  low-frequency  flows.  However,  time-filtering  the  velocities  in 
an  actively  running  model  is  very  expensive  computationally.  Fur¬ 
thermore,  applying  wave  drag  to  only  one  portion  of  the  frequency 
spectrum  of  velocity  can  lead  to  numerical  problems  (see,  e.g.,  Ar¬ 
bic  et  al.,  2010).  Though  the  numerical  problems  can  be  fixed  with 
some  effort  in  topographic  smoothing  (Shriver  et  al.,  2012),  for 
simplicity,  we  assume  that  the  abyssal  flows  in  HYCOM  are  domi¬ 
nated  by  low-frequency  motions.  Some  justification  for  this 
assumption  comes  from  Arbic  et  al.  (2009;  see  their  Section  6), 
where  it  is  shown  that  abyssal  flows  in  eddying  models  that  are 
not  forced  by  tides,  as  is  the  case  in  the  present  manuscript,  are 
dominated  by  low  frequencies  to  a  much  greater  extent  than  are 
abyssal  flows  in  current  meter  observations.  It  is  also  worth  noting 
that  Garner  (2005)  used  his  wave  drag  scheme  on  the  full  velocities 
in  atmospheric  models,  not  on  a  low-passed  component  to  the 
velocities.  In  like  manner,  Nikurashin  and  Ferrari  (2011)  and  Scott 
et  al.  (2011)  did  not  filter  their  ocean  model  velocities  before 
applying  their  wave  drag  schemes  offline.  Thus,  there  is  a  prece¬ 
dent  for  applying  lee  wave  drag  schemes  to  flows  that  include 
some  high-frequency  motions.  Section  3.3  describes  the  C05 
scheme  we  implement  here  for  the  inline  wave  drag  simulation. 
Appendix  A  describes  how  several  parameters  for  the  wave  drag 
parameterization  are  estimated.  Appendix  B  describes  how  wave 
drag  is  extended  to  regions  beyond  those  with  abyssal  hills. 

3.  Bottom  drag  and  wave  drag  schemes 

In  this  section,  we  describe  the  bottom  drag  scheme  used  in  our 
model.  We  then  compare  two  wave  drag  schemes  in  an  offline 
analysis.  The  quadratic  bottom  boundary  layer  drag  and  G05  wave 
drag  schemes  are  both  implemented  inline  in  our  model’s  momen¬ 
tum  equations,  which  will  be  described  in  Section  4.1. 

3.1.  Quadratic  bottom  boundary  layer  drag 

Motivated  by  results  from  idealized  quasi-geostrophic  turbu¬ 
lence  models  on  the  importance  of  either  linear  (Arbic  and  Flierl, 
2004)  or  quadratic  (Arbic  and  Scott,  2008)  bottom  drag  to  the  sta¬ 
tistics  of  mesoscale  eddies,  Sen  et  al.  (2008)  and  Arbic  et  al.  (2009) 
computed  an  energy  dissipation  rate  due  to  bottom  boundary  layer 
drag  of 

EBD  =  pCd\ub\3 .  (1) 

Here,  p  is  the  bottom  density  of  seawater,  Cd  is  the  quadratic  drag 
coefficient  (set  to  0.0025  in  HYCOM),  and  |u6|  is  the  magnitude  of 
the  flow  in  the  bottom  layer.  Although  the  dissipation  rate  is  cubic 
in  the  bottom  current  magnitudes,  the  bottom  stress,  ?BD,  in  the 
momentum  equations  is  given  by, 

?BD  =  -Cd\ub\ub,  (2) 

which  is  quadratic  in  the  velocity  field;  hence  the  convention  of 
referring  to  Cd  as  the  quadratic  drag  coefficient. 
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Fig.  2.  The  log10  of  the  bottom  stratification,  N,  where  N  has  units  [s  '],  at  each  grid  point  from  (a)  the  World  Ocean  Atlas  (WOA)  and  from  (b)  an  average  over  the  thirteenth 
year  of  the  spin-up  period  from  1/12°  HYCOM  without  wave  drag.  Also  shown  is  (c)  the  bathymetic  field,  which  has  units  [m]. 


Quadratic  bottom  boundary  layer  drag  has  a  long  history  in 
oceanography,  dating  back  to  Taylor  (1919).  Evidence  that  the 
parameterization  in  (1)  is  reasonable  comes  from  Egbert  et  al. 
(2004),  who  demonstrate  that  the  dissipation  rates  in  the  shallow 
seas  of  a  forward  global  tide  model  utilizing  quadratic  drag  agree 
quite  well  with  those  inferred  from  satellite  altimetry  data.  Note 
that  (1 )  is  meant  to  represent  the  boundaiy  layer  physics  in  regions 
with  non-negligible  velocities  and  relatively  smooth  topography. 
In  stratified  regions  with  rough  topography,  an  additional  parame¬ 
terization  -  one  for  wave  drag  -  is  required.  Before  delving  into  the 
details  of  the  wave  drag  parameterizations,  we  first  describe  one  of 
their  key  inputs:  a  statistical  representation  of  the  bottom  topogra¬ 
phy  in  regions  where  there  are  abyssal  hills. 

3.2.  Abyssal  hill  rough  topography  representation 

In  order  to  implement  wave  drag,  we  require  a  representation 
of  the  bottom  topography  in  the  spectral  and  physical  domains. 
Goff  and  Jordan  (1988)  demonstrated  the  applicability  of  a  two- 
dimensional  representation  of  the  von  Karman  statistical  model 
for  abyssal  hill  morphology,  the  primaiy  component  of  small-scale 
seafloor  roughness.  The  von  Karman  statistical  model  is  a  band 
limited  fractal  representation,  with  a  power  law  form  at  wavenum¬ 
bers  higher  than  the  corner  wavenumber,  and  flat  below  it.  In  the 
anisotropic  spectral  form  proposed  by  Goff  and  Jordan  (1988),  the 


spectral  representation  of  abyssal  hill  roughness,  P(k,  l ),  is  specified 
by  five  parameters:  root  mean  square  height  ( h );  the  power  law 
exponent  (v,  also  identified  as  the  Hurst  exponent);  corner  wave- 
numbers  in  the  strike  and  normal-to-strike  direction  ( ks  and  k„, 
respectively);  and  the  azimuth  of  the  strike  direction  (£s): 

P(k,  l)  =  47tvh2|Q|-1/2(Yj(k,  0  +  ir(V+1),  (3) 

where 

Q.  =  k2neneTn  +  k2seseTs  (4) 

(with  orthogonal  unit  vectors,  ejjes  =  0),  e„  and  es  are  normal  to  and 
along  the  strike  directions  respectively,  and 

MK 1 )  =  cos^ -&)+(£)  sin2(J-y,  (5) 

where  c  is  the  azimuth  of  the  wavenumber  vector  (fc,/);  hence,  the 
dependence  of  (5)  on  both  components,  k  and  /  (the  zonal  and 
meridional  directions  used  in  Eqs.  (3)-(5)),  of  the  wavenumber  vec¬ 
tor.  Recently,  Goff  and  Arbic  (2010)  estimated  abyssal  hill  statistical 
parameters  globally  using  empirical  relationships  derived  previ¬ 
ously  between  spreading  rate  and  direction,  as  well  as  the  smooth¬ 
ing  effects  of  sediment  cover.  Goff  (2010)  then  formulated  an 
alternative  representation  of  global  abyssal  hill  statistics  based  pri¬ 
marily  on  the  small-scale  roughness  of  the  gravity  field  measured 


124 


D.S.  Trossman  et  al./ Ocean  Modelling  72  (2013)  1 19-142 


by  satellite  altimeters.  Although  both  can  be  considered  realistic 
renderings,  the  latter  is  considered  to  be  more  accurate,  particularly 
in  more  heavily  sedimented  regions  (Goff,  2010).  Goff  (2010)  did 
not,  however,  estimate  azimuthal  orientation.  Thus,  we  will  utilize 
the  parameter  estimates  of  Goff  (2010)  for  h,v,ks,  and  kn,  and  of 
Goff  and  Arbic  (2010)  for  to  generate  the  abyssal  hill  rough 
topography  spectra  used  for  both  the  675  and  G05  wave  drag 
schemes. 

In  theory,  for  the  internal  lee  wave  problem,  we  can  syntheti¬ 
cally  generate  a  topographic  field  of  arbitrary  resolution  using  (3) 
to  arrive  at  a  roughness  input  for  our  wave  drag  estimates.  How¬ 
ever,  the  size  of  such  a  synthetic  map  is  prohibitively  large.  To 
see  this,  we  employ  typical  values  of  the  Coriolis  parameter, 
/~10_4s_1;  ocean  velocity,  LI  —  10  2  to  lms4;  and  abyssal 
Brunt-Vaisala  frequency,  N~10_3s_1  in  a  simple  scaling  argu¬ 
ment  for  the  internal  waves  to  not  be  evanescent, 

f/U  ~  1(T4  itt1  <  \k\  <  N/U  ~  lfT1  m  '.  (6) 

In  (6),  we  used  l/~lms4  for  the  lower  bound  and  U  ~  1CT2  m  s_1 
for  the  upper  bound.  Examination  of  the  highest  wavenumbers 
arising  from  (6)  suggests  that  we  would  need  a  topographic  field 
of  resolution  10  m  or  1/10000°  for  this  task.  In  order  to  keep  the  size 
of  our  topographic  fields  manageable,  we  seek  an  alternative  formu¬ 
lation  (with  computations  done  in  spectral  space)  which  does  not 
require  such  high  resolution  inputs.  Both  the  B75  and  G05  schemes 
are  useful  in  this  sense.  It  should  be  noted,  however,  that//!/  and 
N/U  are  spatially  varying  quantities.  We  cannot  know  what  the 
stratification  or  velocity  field  will  be  at  each  horizontal  grid  point 
as  a  function  of  time  before  we  perform  our  model  simulations. 
Thus,  we  need  to  impose  an  assumption  on  the  range  of  relevant 
wavenumbers  via  (6). 

3.3.  Wave  drag  schemes 


over  a  topographic  feature)  in  an  ad  hoc  way  via  (9).  This  motivates 
an  alternative  wave  drag  scheme  that  we  describe  below,  and  that 
has  been  successfully  used  in  an  altered  form  (assuming  oscillating 
background  flows)  in  a  prognostic  ocean  tide  model  (Arbic  et  al., 
2004). 

G05  formulated  a  scheme  to  parameterize  the  sink  of  horizontal 
momentum  due  to  the  interaction  of  a  stratified  atmosphere  with 
orography.  Using  the  same  scheme,  we  can  parameterize  the 
momentum  sink  that  occurs  when  low-frequency  oceanic  flows 
impinge  upon  bottom  topography.  Adding  this  term  to  the 
momentum  equations  will  lead  to  a  sink  in  the  energy  as  well 
(see  discussion  below).  In  contrast  to  the  B75  scheme,  the  G05 
scheme  is  based  on  a  calculation  of  the  column-integrated  momen¬ 
tum  forcing  that  is  exact  in  the  small-amplitude  (linear)  limit. 
Since  orography  is  never  entirely  linear,  the  G05  scheme  adjusts 
the  forcing  for  partially  blocked  or  partially  deflected  flow  based 
on  a  dimensional  analysis.  An  important  assumption  is  that  the 
blocked  part  of  the  flow  adjusts  independently  for  each  topo¬ 
graphic  feature.  This  is  expected  to  be  less  accurate  for  more  den¬ 
sely  populated  regions  of  topographic  features.5  The  scheme  uses 
additional  assumptions  to  distribute  the  forcing  in  the  vertical  col¬ 
umn.  For  this  preliminary  study,  we  replace  these  assumptions  with 
a  simpler  one  that  all  momentum  flux  is  deposited  in  the  lowest 
500  m.  The  implications  of  a  forcing  that  extends  to  higher  levels 
are  left  for  a  future  study. 

To  calculate  the  base  flux  of  momentum  according  to  the  G05 
scheme,  we  need  the  information  tensor  T  (units  [kg  m  2  s  1  ]), 
describing  the  specified  subgrid  topography.  The  information  ten¬ 
sor  can  be  calculated  from 


T(x,y) 


■^r/dM|P(k,0lg 

IW)lf 


^fdkdl\P(k,l)\ 
£?fdkdl \P(k,l)\ 


Tr,! 

Tl  .2 

T2,  i 

T2, 2. 

(10) 


B75  formulated  a  scheme  that  relates  a  spectral  representation 
of  the  topographic  field’s  abyssal  hill  roughness,  P(k,l),  to  the  en¬ 
ergy  flux  into  lee  waves  from  the  geostrophic  flow,  U,  impinging 
upon  bottom  topography.  The  energy  dissipation  rate  per  unit  area 
can  be  expressed  as 

Eveii  =  ^  I"'1"'  dkP(k) \J (N2  -  \U\2k2)(\U\2k2 -P),  (7) 

27t  J |/|/|U| 

(675)  where 

P(k)=-L  r  dlBp(k,l),  (8) 

2 n  J-0 „  |k| 

for  k  =  ( k ,  /),  the  wavenumber  vector  in  the  reference  frame  in 
which  k  is  along  and  /  is  across  the  mean  flow,  U,  respectively. 
The  range  of  wavenumbers  used  in  (7)  is  determined  by  (6).  This 
wave  drag  scheme  can  be  used  offline  on  model  output,  as  in  Nik- 
urashin  and  Ferrari  (2011)  and  Scott  et  al.  (2011).  The  latter  study 
made  use  of  a  multiplicative  correction  factor,  Efac,  to  the  675 
scheme’s  energy  dissipation  rate,  (7).  The  correction  factor  is  a  func¬ 
tion  of  the  Froude  number  of  the  bottom  flow,  Fr  =  U/(NZ),  where  Z 
is  a  vertical  length  scale: 

Efac  =  L  [cos-'(l  -  21)  -  2(1  -  21)^1  -L)]  ■  (9) 

Here,  L  =  1  when  Fr  '  ^  0.751  and  L  =  0.75 Fr  otherwise.  We  fol¬ 
low  suit  in  our  offline  675  computations,  with  the  intent  to  compare 
our  results  with  estimates  from  Scott  et  al.  (201 1 ).  The  675  scheme, 
and  its  implementation  by  Scott  et  al.  (2011),  relies  on  a  two- 
dimensional  field  that  represents  statistical  features  of  the  underly¬ 
ing  terrain  at  each  grid  point.  Furthermore,  Scott  et  al.  (2011)  only 
accounts  for  topographic  blocking  (i.e.,  flow  that  may  not  propagate 


where  p  and  N  are  the  instantaneous  density  and  Brunt-Vaisala  fre¬ 
quency,  respectively,  each  averaged  over  the  bottom  500  m  (see  be¬ 
low).  We  use  the  local  power  spectrum,  P(k,l),  obtained  from  (3) 
using  the  five  parameters  estimated  from  Goff  and  Arbic  (2010) 
and  Goff  (2010),  to  characterize  subgrid  features  of  the  oceans  abys¬ 
sal  hills  centered  on  the  grid  cell. 

Based  on  the  dispersion  relation  for  internal  waves,  we  expect 
the  scales  relevant  to  the  wave  drag  to  fall  in  the  range  given  by 
(6),  which  are  used  to  compute  T  via  (10).  Scott  et  al.  (2011),  by 
contrast,  used  the  range,  1CT6  <  |k|  <  1CT1  m  '.  However,  wave- 
numbers  between  10_6m_1  and  10_4m_1  may  only  be  relevant 
close  to  the  equator  or  at  mid-latitudes  where  either/ is  very  small 
or  U  is  relatively  large,  making  the  left-hand  side  of  (6)  small.  Using 
the  same  wavenumber  range  as  Scott  et  al.  (201 1 )  at  most  doubles 
the  magnitudes  of  the  components  of  T  in  a  few  scattered  locations 
such  as  along  the  Mid-Atlantic  Ridge.  Because  of  feedbacks  be¬ 
tween  velocities  and  wave  drag,  doubling  T  would  result  in  less 
than  a  factor  of  two  increase  in  the  dissipation  (see,  for  instance, 
the  ~log10(2)  reduction  with  a  factor  of  two  increase  in  drag  in 
Fig.  4a  of  Arbic  et  al.  (2004)  using  a  forward  one-layer  ocean  tide 
simulation). 

As  in  G05,  we  assume  that  the  total  drag  is  in  the  direction  of  the 
linear  drag,  ?lin  =  T ud,  and  that  its  magnitude  is  determined,  via  the 
aforementioned  dimensional  scaling,  by 


5  Another  important  assumption  is  that  the  G05  scheme  does  not  explicitly 

consider  the  effects  of  rotation  on  the  drag  itself,  the  assumption  being  that  N  >/. 
However,  the  information  tensor  used  to  calculate  the  wave  drag  is  integrated  over 
the  range  of  wavenumbers  specified  by  (6),  which  does  depend  upon  the  Coriolis 
parameter. 
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?  =  (Xx,  Xy) 


{w+w) 


(Jud). 


(11) 


Here,  ud  is  the  velocity  averaged  over  the  bottom  Hw D  =  500  m,  Dp 
and  D„p  (respectively)  are  the  propagating  and  non-propagating 
parts  of  the  drag  based  on  the  dimensional  analysis,  and  D*  is  the 
linear  limit  of  Dp.  Our  choice  of  HWD  =  500  m  is  guided  by  the  ver¬ 
tical  decay  scale  in  St.  Laurent  et  al.  (2002),  which  is  based  on  the 
microstructure  observations  of  St.  Laurent  et  al.  (2001)  in  the  Brazil 
Basin.  A  more  recent  observational  study  (St.  Laurent  et  al.,  2012) 
was  done  in  the  Drake  Passage  and  found  no  evidence  of  enhanced 
turbulence  at  heights  exceeding  500  m  above  the  bottom.  As  in  HY- 
COM  simulations  with  embedded  tides  (Arbic  et  al.,  2010), 
Hwd  =  500  m  is  much  greater  than  HBD  =  10  m,  the  thickness  over 
which  quadratic  bottom  drag  acts. 

The  expressions  for  Dp,  Dnp,  and  D*  are  taken  from  G05  without 
modification  -  his  Eqs.  15,16.  For  the  convenience  of  the  reader, 
they  are  given  here: 


D*  =  a0 


P^Hy \(.2y-t)(H2m+r~H2+r) 

NLr  r  (2  +y-e)(H2^-H2Zne) 


-n.wPVi 


Dp  =  a0H\ 


NLr  H2y-e  -  H2y ~ 

'  1  1  mnv  1  ■*  min 


Tj2+y-e  _  L/2+y- 
flclip  flmin 


Dnp  —  Q\H'r 


2  +  y-e 
pV3 


ur-e-fi  _  ur-e-P 

I_f2+p  max  “clip 


y-e~, 


2y-e 


NLr(  1  +  p)  H2y-e  -  H2y~e 

'  v  1/  ‘  lmnx  1 1  mm 


Tii+y-e  _  ul+y-e  uy-e-P  _  uy-e-P' 

nrnax  nclip  +pnmax  nclip 

1 +y-e  crit  y-e-p 


(12) 


The  minimum  and  maximum  terrain  heights,  Hmin  =  0AHmax  and 
Hmax,  are  computed  in  the  way  suggested  by  G05  and  normalized 
by  V/N,  the  internal  scale.  For  this  scale,  we  use  the  velocity,  V, 
and  stratification,  N2,  in  the  lowest  HWD  =  500  m  of  the  model. 
The  terrain  height  limit  that  prevents  wave  drag  from  increasing  be¬ 
yond  a  certain  extent  is  Hdip  =  min{Hmax,  max{Hmin. Hcrjt}},  all  vari¬ 
ables  with  the  subscript  ‘r’  (for  reference)  cancel  out,  and 
Kent  =  NHait/V  is  the  critical  nondimensional  height  that  deter¬ 
mines  the  degree  to  which  the  flow  is  blocked.  The  critical  nondi¬ 
mensional  height  is  a  universal  parameter  in  this  scheme,  leaving 
Hcrit  to  be  flow-dependent.  We  set  Hcrit  =  0.7.  Other  choices  of 
parameters  used  in  the  G05  scheme,  such  as  the  coefficients  for 
the  propagating  and  non-propagating  components  of  wave  drag 
(a0  and  a, ,  respectively),  y,  e,  and  p  are  discussed  in  Appendix  A. 

For  simplicity,  and  for  a  more  direct  comparison  with  the  B75 
scheme,  we  reduce  the  parameterized  tensor  drag  in  the  G05 
scheme  to  a  scalar  drag  with  the  same  instantaneous  impact  on 
the  kinetic  energy.  This  procedure  has  been  used  in  tidal  simula¬ 
tions  (e.g.,  Arbic  et  al.,  2010).  For  the  stress,  x  ,  given  in  (11)  we 
compute 


Xdrag  — 


X  -Ud 
P\Ud\2’ 


(13) 


which  is  a  decay  rate  times  a  vertical  length  scale  (units  [ms-1]). 
The  value  of  rdrag  is  computed  using  the  full  tensor,  T,  at  each  grid 
point  inline  with  HYCOM.  To  be  precise,  the  stress  associated  with 
the  wave  drag  term  in  the  momentum  equations  is  computed  as 


XwD  =  -|tdrng|Ud- 


(14) 


Table  1 

The  estimated  globally  integrated  energy  dissipation  rate  due  to  wave  drag 
[TW  =  1012  W]  using  two  different  wave  drag  schemes:  B75  (Bell)  and  G05  (Garner). 
We  utilize  a  five-day  averaged  velocity  field,  averaged  over  the  bottom  500  m,  from 
the  final  year  of  the  1  /12°  HYCOM  spin-up  without  wave  drag.  Two  different  bottom 
stratification  fields  are  used:  from  the  World  Ocean  Atlas  (WOA)  and  an  average  over 
the  bottom  500  m  from  the  final  year  of  the  1/12°  HYCOM  spin-up  without  wave 
drag. 


Wave  drag  scheme 

WOA 

HYCOM 

Garner 

0.47 

0.45 

Bell 

0.45 

0.51 

the  velocities  and  stratification  in  the  model’s  bottom  500  m  in 
addition  to  features  of  the  underlying  topography. 

The  use  of  (14)  slightly  changes  the  time-dependent  flow  com¬ 
pared  to  the  full  tensor  implementation  in  Arbic  et  al.  (2004)  and 
G05,  in  which  the  drag  is  not  constrained  to  be  parallel  to  the  flow. 
The  angle  between  the  flow  and  drag  is  mainly  determined  by  the 
difference  between  the  diagonal  elements  of  T  (the  off-diagonal 
elements  are  inherently  small).  In  Section  4.2,  we  will  explicitly 
compare  x/p  with  rdmgud. 

4.  Energy  budgets 

In  this  section  we  list  the  source  and  sink  terms  in  the  total 
mechanical  energy  budget.  We  follow  this  with  a  comparison  of 
wave  drag  contributions  to  the  energy  budget  using  the  B75  and 
G05  schemes  acting  offline  on  output  from  HYCOM  simulations 
that  do  not  have  wave  drag  actively  altering  the  velocities.  We  also 
compare  offline  estimates  using  bottom  stratification  fields  from 
both  WOA  (World  Ocean  Atlas  2005)6  and  HYCOM,  the  latter  calcu¬ 
lated  with  the  TEOS-10  package  (McDougall  and  Barker,  2011).  We 
compute  the  input  and  output  terms  of  an  energy  budget  in  line 
for  HYCOM  simulations  with  and  without  wave  drag,  the  former 
including  the  G05  wave  drag  scheme  actively  altering  both  the  bot¬ 
tom  velocities  and  bottom  stratification  fields. 

4.1.  Total  mechanical  energy  budget 

The  kinetic  energy  equation,  part  of  the  total  mechanical  energy 
budget,  here  is  derived  from  the  momentum  equations, 

^  +  (u  ■  V)uH  +  - Vp  +fk  x  uH  +  gk 
at  p 

=  9<"Z  >  ~Hs'>  -  0{z  <  Hbd  -  zb)  |uH|uH  -  B{z 

P  ns  tiBD 

<  HWD  -  Zb)  ^  UH  -  4;  ( Vz  4  fa)  ~  V  •  ( Vh,2  VtiH 

Hwd  oz  \  oz  j 

+  vmVV2Sh),  (15) 

where  V  =  (d/dx,d/dy,d/dz)  is  a  three-dimensional  gradient  oper¬ 
ator,  uH  =  (u.v)  is  the  two-dimensional  velocity  along  isopycnal 
surfaces,  fi  =  (u,  v ,  w)  is  the  three-dimensional  velocity,  p  is  the 
pressure,  k  is  a  unit  vector  in  the  vertical  direction,  g  =  9.806  m  s~2 
is  the  acceleration  due  to  gravity,  and  /  is  the  Coriolis  parameter. 
The  right  hand  side  of  (15)  displays  the  wind  stress,  bottom  drag, 
wave  drag,  vertical  viscosity,  and  horizontal  viscosity  terms,  respec¬ 
tively.  6(z  <  HBd  -  Zb)  is  a  step  function  that  is  one  in  the  layer  of 
thickness  HBD  =  10  m  from  the  bottom,  and  zero  for  all  other  layers. 
0(z  <  HWd  -  zb)  is  a  step  function  that  is  one  in  the  layer  of  thickness 
HWd  =  500  m  from  the  bottom,  and  zero  for  all  other  layers. 
9(z  >  -Hs)  is  a  step  function  that  equals  one  in  the  surface  layer  of 


Thus,  in  the  model  momentum  equations,  |rdrag|  in  (14)  is  to  wave 
drag  as  C^Uj,!  in  (2)  is  to  bottom  drag.  While  Cdjiii,|  depends  only 
on  the  velocities  in  the  model’s  bottom  layer,  \rdmg\  depends  on  both 


6  http://www.nodc.noaa.gov/OC5/WOA05/woa05data.html 

7  In  regions  shallower  than  500  m,  Hw0  is  taken  to  be  the  depth  of  the  water 
column. 
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thickness  Hs  and  zero  for  all  other  layers.  twM  is  the  surface  wind 
stress,  which  is  computed  accounting  for  the  relative  velocities  of 
the  wind  and  surface  currents  (Scott  and  Xu,  2009).  The  vertical  vis¬ 
cosity,  vz  ~  1CT4  -  1CT3  m2  s-1,  is  calculated  inline  in  HYCOM  using 
KPP  (Large  et  al.,  1994).  The  horizontal  viscosity,  which  is  on  the  or¬ 
der  of  102  -  103  m2  s  \  is  also  calculated  inline  in  HYCOM  with  vh2 
taken  to  be  the  maximum  of  a  Laplacian  term  and  a  Smagorinsky 
(1993)  term  and  vh4  taken  to  be  a  biharmonic  term. 

Multiplying  (15)  by  density,  dotting  that  equation  with  the 
velocity,  and  performing  a  volume  integral,  one  finds  the  kinetic 
energy  equation,  which  includes  a  conversion  term  between  ki¬ 
netic  energy  and  potential  energy.  We  first  write  the  kinetic  energy 
equation  here  for  an  incompressible,  hydrostatic  fluid  explicitly  as 


P EKtime  T"  PppUthi  —  P pressure  P  input  P output  “F  QlK->Ep 


(16) 


(Griffies,  2004)  where  the  area  integrals  are  evaluated  at  the 
surface  and  arise  as  a  result  of  invoking  Gauss’  Theorem, 
Ek  =  pu-  u/2  is  the  kinetic  energy  per  unit  volume,  Pinput  =  Pwmd, 
P output  =  Pbd  +  Pwd  +  Pw  +  Phv .  and 
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J  dx  J  dy(p|rdmg||ud|2 


Pbd  — 

Pwd  = 

Pw  =  I  dx 


Phv  = 


rl 


-  d  .  d  _  , 

pU¥z^dzUH) 


(17) 


dy{L 

j  dx  I  dy(^  J  dz^pu-  V(vhj2VuH  +  vMVV2uH)] 
*->e,  =  -  J  dx  J  dy  (/z  dz[gpw\  j  . 


PEKtime  denotes  the  time  rate  of  change  in  globally  integrated  kinetic 
energy,  PEKadi>  denotes  the  change  in  kinetic  energy  due  to  advective 
fluxes  across  the  surface,  Ppressure  denotes  the  divergence  of  kinetic 
energy  due  to  pressure  differentials  at  the  surface  (evaluated  at 
the  atmospheric  pressure  p  =  patm  which  is  constant  in  the  model), 
Cek->e p  denotes  the  conversion  of  kinetic  energy  ( EK )  to  potential 
energy  (£P),  Pwinc /  denotes  the  wind  power  input  to  the  total  flow, 
Pbd  denotes  the  kinetic  energy  dissipation  due  to  bottom  drag, 
Pwd  denotes  the  kinetic  energy  dissipation  due  to  wave  drag,  Pw  de¬ 
notes  the  kinetic  energy  dissipation  due  to  vertical  viscosity,  and 
Phv  denotes  the  kinetic  energy  dissipation  due  to  horizontal  viscos¬ 
ity,  respectively.  Here,  t;  is  the  sea  surface  height,  f  dxf  dy(-)  indi¬ 
cates  an  integral  over  the  surface  area  of  the  globe  and  at  a  given 
horizontal  grid  point  in  HYCOM,  zb  is  the  depth  of  the  ocean,  ub  is 
the  two-dimensional  velocity  averaged  over  the  bottom  HBD  meters, 
ud  is  the  two-dimensional  velocity  averaged  over  the  bottom  Hwo 
meters,  us  is  the  two-dimensional  surface  velocity,  and  w  is  the  ver¬ 
tical  velocity.  In  addition  to  the  wind  power  input  to  the  total  flow, 
we  compute  the  wind  power  input  to  the  geostrophic  flow  via 
f  dxf  dy(Ug  ■  TWjnd),  where  ug  is  the  geostrophic  velocity  derived 
from  the  sea  surface  heights  (Wunsch,  1998),  in  order  to  evaluate 
how  wave  drag  impacts  the  surface  in  Section  4.3. 

It  is  worth  noting  that  the  integrands  of  both  the  horizontal  and 
vertical  viscosity  terms  in  (17)  can  represent  energy  sources  (i.e., 


have  positive  values)  at  isolated  grid  points  (see  Section  4.3).  For 
horizontal  viscosity,  this  may  seem  to  deviate  from  the  spectral 
model-based  argument  that  the  Laplacian  term  is  a  net  dissipator. 
However,  the  latter  argument  only  applies  to  an  integral  over  a 
domain  much  larger  than  a  single  grid  point.  Fox-Kemper  and 
Menemenlis  (2008)  note  that  the  Smagorinsky  (1993)  parameteri¬ 
zation  takes  a  domain-averaged  value  of  the  horizontal  viscosity 
and  replaces  it  with  a  local  value.  Therefore,  the  small-scale  spatial 
variation  in  maps  of  energy  dissipation  by  the  horizontal  viscosity 
term  should  not  be  interpreted  as  physical.  The  small-scale  spatial 
variation  of  the  vertical  viscosity  term  may  not  have  a  physical 
interpretation  either.  The  global  integrals  and  some  regional  char¬ 
acteristics  of  the  horizontal  and  vertical  viscosity  terms  in  (17), 
however,  are  physically  meaningful. 

The  potential  energy  equation,  also  part  of  the  total  mechanical 
energy  budget,  is  derived  by  making  use  of  the  equation  (Griffies, 
2004),  dp/dt  =  V  ■  (tcVp),  where  k  is  the  eddy  diffusivity  tensor; 
the  definition  of  the  vertical  velocity,  dz/dt  =  w;  the  potential  en¬ 
ergy  per  unit  volume,  EP  =  pgz\  the  definition  of  a  material 
derivative, 


fdxfdyfdz'ff.ftofdyf 
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and  Gauss’  Theorem, 
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Going  from  left  to  right,  the  terms  on  the  right-hand  side  of  (18)  are 
the  time  rate  of  change  in  potential  energy  and  the  reversible 
changes  in  potential  energy  via  advective  fluxes  at  the  surface.  Sim¬ 
ilarly  going  from  left  to  right,  the  terms  on  the  right-hand  side  (last 
line)  of  (19)  are  the  conversion  between  potential  energy  and  ki¬ 
netic  energy  via  buoyancy  fluxes,  the  rate  of  change  of  potential  en¬ 
ergy  due  to  diffusive  mass  flux  that  acts  to  mix  or  stir  density 
(evaluated  at  the  surface  and  at  the  seafloor)  across  isopycnals, 
the  rate  of  change  of  potential  energy  due  to  diffusive  mass  flux  that 
acts  to  mix  or  stir  density  along  isopycnals,  and  the  conversion  from 
internal  energy  to  potential  energy.  We  assume  that  the  diffusive 
mass  flux  (i.e.,  buoyancy  diffusion)  evaluated  at  the  seafloor  is  neg¬ 
ligible  because  there  is  no  geothermal  heat  flux  there  in  the  model. 
We  also  neglect  the  along-isopycnal  contributions  to  the  buoyancy 
diffusion  term  by  assuming  that  VHp  =  (d/dx,d/dy)p  =  0  in  the 
model,  reducing  K  to  a  scalar.3  We  allow  for  the  diffusivity,  k  =  vz, 
at  the  surface  to  spatially  vary  in  the  diffusive  mass  flux  across  iso¬ 
pycnals  and  conversion  from  internal  energy  to  potential  energy 
terms  on  the  right-hand  side  of  (19),  but  (19)  is  otherwise  the  same 
as  Eq.  (6)  in  Winters  et  al.  (1995).  The  vertical  velocity  at  the  surface, 
needed  for  PEpadv  and  Pnguswet  is  calculated  by  making  use  of  the  sur¬ 
face  boundary  condition  that  w  =  drj/dt  +  uH  ■  Vi).  We  do  not  in¬ 
clude  work  done  associated  with  the  compressibility  of  fluid  here 
since,  by  our  assumptions,  V  ■  u  =  0.  The  buoyancy  flux,  or  thermo¬ 
dynamic  work,  term  at  the  surface  does  not  explicitly  appear  in  the 
conservation  equation  because,  as  argued  by  Wunsch  and  Ferrari 


8  When  VhP  ^  0,  it  is  unclear  how  to  compute  the  along-isopycnal  diffusivities  as 
HYCOM  runs  because  the  Laplacian  and  biharmonic  contributions  to  the  horizontal 
viscosity  have  diffusion  velocities  from  which  we  can  compute  diffusivities,  but  the 
Smagorinsky  term  does  not. 
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(2004),  the  buoyancy  flux  at  the  surface  enters  as  internal  energy 
which  then  can  be  converted  to  other  forms  of  energy  by  changing 
the  mass  of  the  ocean  via  the  last  term  in  (19).  Saenz  et  al.  (2012) 
suggest  that  the  available  potential  energy  ultimately  provided  by 
the  buoyancy  flux  at  the  surface  can  be  converted  into  kinetic  energy 
via  the  first  term  on  the  right-hand  side  of  (19).  Tailleux  (2010)  sug¬ 
gests  that  the  buoyancy  flux  term  may  be  a  significant  contributor  to 
the  energy  budget,  and  von  Storch  et  al.  (2012)  found,  using  their 
high-resolution  ocean  model,  that  the  contribution  of  the  buoyancy 
flux  to  the  input  of  potential  energy  is  comparable  to  the  wind  power 
input. 

By  equating  (18)  with  (19),  we  write  the  potential  energy  equa¬ 
tion  succinctly  here  as 

P Eptime  “F  P EPadv  ~  P diffusive  '  f-£p->£K  “F  (-Ei  :>Ep ,  (20) 

where 


'  Eptime  — 


=  dx 


f  r  dz 

dEp 

[tftj 

PEpadv  =  I  dx  J  dy(ijgpw) 

'diffusive  =  J  dx  j  dy(gi]K^ 

:p->Ep  =  J  dX  f  dy[J_z  dZ{gPW 1 


P  (In 

Ce,  . 

Ce,->ep  =  -  /  dx 


HI 


dz 


gK 


(21) 


PEptime  denotes  the  time  rate  of  change  in  globally  integrated  poten¬ 
tial  energy,  PEpadv  denotes  the  changes  in  potential  energy  due  to 
advective  fluxes  at  the  surface,  Pdiffusive  denotes  the  potential  energy 
change  due  to  diffusive  mass  fluxes  (buoyancy  diffusion),  C£p_>£|r 
denotes  the  conversion  of  potential  energy  to  kinetic  energy,  and 
CEi->ep  denotes  the  conversion  from  internal  energy  (£/)  to  potential 
energy.  The  advective  flux  of  potential  energy  should  be  thought  of 
as  the  work  the  ocean  is  doing  on  the  atmosphere.  Again,  note  that 
here  we  evaluate  the  diapycnal  component  of  the  diffusive  mass 
fluxes.  The  buoyancy  diffusion  term,  Pdiffusive,  can  be  interpreted  as 
the  power  associated  with  density  mixing  or  stirring.  Finally,  our  to¬ 
tal  mechanical  energy  budget  is  calculated  by  summing  (16)  and 
(20),  yielding 


1  EKtime  "F  ‘  Eptime  <  *  tKaav  i  *  tPaav 

—  P pressure  Pdiffusive  "F  P- EI->EP  "F  P input 


(22) 


Each  of  the  advective  flux  terms  and  the  pressure  term  make  use  of 
the  incompressibility  assumption. 

It  is  important  to  note  that,  with  the  exception  of  the  partial 
time  derivative  terms,  the  terms  in  (22)  in  our  inline  calculations 
are  evaluated  instantaneously  at  each  baroclinic  time  step  of  two 
minutes.  For  example,  the  average  of  us  ■  x^d  is  computed  instead 
of  the  product  of  the  averages  of  us  and  x^d.  Thus,  we  are  confi¬ 
dent  that  aliasing  higher  frequency  motions  forced  by  2003  NO¬ 
GAPS  winds  is  not  a  problem.  The  partial  time  derivative  terms, 
P Eptime  and  PEptime,  are  diagnosed  using  the  global  integrals  of 
monthly  averaged  fields  from  the  final  year  of  the  HYCOM  simula¬ 
tion  with  wave  drag. 


4.2.  Preliminary  offline  estimates  of  wave  drag 

In  this  section,  we  perform  offline  estimates  of  energy  dissipa¬ 
tion  rates  due  to  wave  drag.  First,  we  show  that  the  G05  scheme 
yields  spatial  maps  and  globally  integrated  dissipation  rates  that 
are  comparable  to  those  computed  from  the  B75  scheme.  To  do 
this,  we  use  five-day  averaged  velocities,  averaged  over  the  bottom 


500  m,  from  the  fourteenth  year  of  a  1/12°  HYCOM  simulation 
spun-up  without  wave  drag.  We  use  two  fields  for  the  mean  annual 
abyssal  stratification:  one  taken  from  the  same  model  simulation, 
and  the  other  taken  from  WOA.  The  WOA  stratification  (Fig.  2a) 
is  typically  weaker  than  the  HYCOM  stratification  (Fig.  2b),  partic¬ 
ularly  in  areas  of  rough  abyssal  hills  that  flank  topographic  ridges 
(Fig.  2c).  For  the  offline  estimates,  we  also  use  the  abyssal  hill 
rough  topography  power  spectra  that  were  computed  with  param¬ 
eters  from  Goff  (2010)  and  Goff  and  Arbic  (2010). 

The  offline  globally  integrated  energy  dissipation  rates  due  to 
wave  drag  are  listed  in  Table  1,  while  spatial  maps  are  shown  in 
Fig.  3.  When  the  WOA  bottom  stratification  field  is  used,  the  glob¬ 
ally  integrated  energy  dissipation  rate  using  the  C05  scheme  is  lar¬ 
ger  than  that  using  the  B75  scheme  by  about  5%.  On  the  other  hand, 
the  globally  integrated  energy  dissipation  rate  using  the  C05 
scheme  is  smaller  than  those  using  the  B75  scheme  by  about  10% 
when  the  HYCOM  stratification  field,  averaged  over  the  bottom 
500  m,  is  used.  These  differences  are  likely  because  the  B75  scheme 
does  not  explicitly  distinguish  between  flow  that  propagates  and 
flow  that  does  not  propagate  over  a  topographic  feature,  whereas 
the  G05  scheme  does.  The  nonlinearity  of  the  G05  scheme  depends 
on  the  stratification,  and  the  B75  scheme  depends  differently  upon 
the  stratification.  It  is  also  possible  that  the  differences  are  due  to 
the  explicit  dependence  of  the  B75  estimates  upon  the  Coriolis 
parameter.  However,  the  discrepancies  do  not  appear  to  be  related 
to  the  Coriolis  parameter  because  differences  between  the  upper 
panels  (G05  scheme)  and  lower  panels  (B75  scheme)  in  Fig.  3  do 
not  appear  to  be  Iatitudinally  dependent  any  more  than  they  are 
longitudinally  dependent. 

The  regions  with  the  greatest  energy  dissipation  rates  using 
either  wave  drag  scheme  are  in  the  Subpolar  North  and  South 
Atlantic  Oceans,  the  western  Indian  Ocean,  and  the  Southern  Ocean 
(Fig.  3).  The  B75  scheme  has  larger  dissipation  rates  than  the  G05 
scheme  off  the  coast  of  Brazil,  in  the  South  Pacific  Ocean,  and  in 
the  western  Indian  Ocean.  On  the  other  hand,  the  B75  scheme 
has  smaller  dissipation  rates  than  the  G05  scheme  in  the  wake  of 
the  Drake  Passage,  in  the  Antarctic  Circumpolar  Current  (ACC),  in 
the  Aghulas  Current  region,  and  in  the  Gulf  Stream  Extension  re¬ 
gion.  Both  the  B75  and  G05  schemes  estimate  greater  dissipation 
rates  (Fig.  3)  in  regions  where  we  would  expect  more  energy  dissi¬ 
pation  due  to  a  combination  of  relatively  strong  flow  (not  shown) 
and  relatively  strong  stratification  (Fig.  2b)  over  rough  topographic 
features  (Fig.  2c). 

Our  offline  dissipation  rates  using  the  WOA  stratification  are 
larger  than  those  of  Nikurashin  and  Ferrari  (2011)  (0.2  TW)  and 
are  closer  to  those  of  Scott  et  al.  (2011)  (0.34-0.49  TW).  We  have 
used  a  similar  WOA  stratification  field  and  the  same  abyssal  hill 
rough  topography  power  spectra  as  the  larger  of  the  estimates  of 
Scott  et  al.  (2011),  but  our  abyssal  velocity  field  is  different.  The 
abyssal  stratification  field,  abyssal  velocity  field,  and  abyssal  hill 
rough  topography  power  spectrum  formulation  are  all  different 
here  from  those  used  in  Nikurashin  and  Ferrari  (2011).  Nikurashin 
and  Ferrari  (2011)  used  an  isotropic  spectral  form  of  the  abyssal 
hills,  while  Scott  et  al.  (2011)  allowed  for  anisotropy  to  enter  into 
the  abyssal  hill  power  spectra.  The  largest  difference  between  the 
estimates  of  Nikurashin  and  Ferrari  (2011)  versus  those  in  Scott 
et  al.  (2011)  and  in  the  present  study  is  along  the  equator,  where 
the  former  study  found  much  larger  dissipation  rates.  Because  each 
of  the  aforementioned  studies  used  similar  stratification  products, 
it  is  likely  that  the  relatively  large  dissipation  rates  found  along  the 
equator  in  the  estimates  of  Nikurashin  and  Ferrari  (2011)  are  due 
to  their  use  of  an  isotropic  spectral  form  of  the  abyssal  hill  power 
spectrum.  However,  the  velocity  field  used  by  Nikurashin  and 
Ferrari  (2011)  could  be  responsible  for  their  relatively  large 
dissipation  rates  found  along  the  equator.  The  stratification  field 
may  also  be  important  for  a  given  velocity  field,  evidenced  by  the 
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(a)  (Gamer  scheme,  N  from  WOA)  (b)  (Garner  scheme,  N  from  HYCOM) 


Fig.  3.  Spatial  maps  of  the  logi0  of  the  offline  estimates  of  energy  dissipation  rates  due  to  wave  drag,  where  the  energy  dissipation  rates  have  units  [W  nr2],  using  two 
different  estimation  schemes:  (a,  b)  GO 5  (Garner)  and  (c,  d)  B75  (Bell),  and  using  two  different  bottom  stratification  fields:  from  the  World  Ocean  Atlas  (WOA)  (a  and  c)  and 
from  an  average  over  the  thirteenth  year  of  the  spin-up  period  from  1  /12”  HYCOM  without  wave  drag  (HYCOM)  (b  and  d).  The  white  regions  are  non-abyssal  hill  regions. 
Listed  in  Table  1  are  the  globally  integrated  energy  dissipation  rates  of  each  estimate  [TW  =  1012  W], 


relatively  large  wave  drag  dissipation  rates  found  here  along  the 
equator  with  the  HYCOM  stratification  but  not  with  the  WOA  strat¬ 
ification  (Fig.  3).  Therefore,  we  will  investigate  whether  differences 
between  velocity  and  stratification  products  can  play  an  important 
role  in  the  offline  wave  drag  dissipation  rate  estimates  along  the 
equator  in  Section  4.4. 

It  is  important  to  note  that,  as  in  Nikurashin  and  Ferrari  (2011) 
and  Scott  et  al.  (2011),  our  offline  estimates  of  wave  drag  are  only 
active  in  regions  where  there  are  abyssal  hills.  A  globally  inte¬ 
grated  offline  estimate  using  the  C05  scheme  (with  either  abyssal 
stratification)  together  with  drag  fields  that  are  also  defined  in 
non-abyssal  hill  regions  yields  about  three  times  as  much  energy 
dissipation  as  estimates  made  only  in  abyssal  hill  regions.  This  is 
almost  entirely  owing  to  large  energy  dissipation  rates  found  in 
the  Drake  Passage,  around  Iceland,  and  in  the  Bering  Strait.  All  of 
these  regions  are  marked  by  either  relatively  strong  flow  or  rela¬ 
tively  strong  stratification  over  rough  topographic  features. 

Before  moving  onto  the  inline  estimates  of  energy  dissipation 
rates  due  to  wave  drag,  we  compare  x/p  from  (11)  with  Ir^gldd 
from  (14)  in  Fig.  4  using  the  same  abyssal  velocity  and  stratifica¬ 
tion  fields  as  in  Fig.  3a.  The  largest  differences  between  x/p  and 
\rdrag\tid  are  not  localized  to  steep  topographic  ridges,  and  these  dif¬ 
ferences  vary  between  different  choices  of  five-day  averaged  abys¬ 
sal  velocity  fields.  However,  Fig.  4  shows  many  of  the  typical 
discrepancies.  The  differences  can  be  particularly  pronounced  in 


the  South  Indian  Ocean  (Fig.  4a  and  b),  east  of  Drake  Passage 
(Fig.  4a  and  b),  around  the  Gulf  Stream  (Fig.  4e  and  f),  and  between 
Great  Britain  and  Iceland  (Fig.  4e  and  f).  These  are  regions  where 
the  off-diagonal  components  of  T  can  be  significant  enough  to  ro¬ 
tate  the  wave  drag  several  degrees  from  the  direction  of  the  abys¬ 
sal  velocities  and  contract  the  length  of  the  wave  drag  vector  by 
several  percent.  The  differences  are  relatively  small  elsewhere, 
such  as  in  the  equatorial  Pacific  Ocean  (Fig.  4c  and  d).  While  it  is 
important  to  note  that  we  have  only  incorporated  the  contribution 
of  wave  drag  to  the  momentum  equations  in  the  direction  parallel 
to  the  flow,  the  inline  calculation  adjusts  the  velocities  to  the  wave 
drag.  The  mechanical  energy  budget,  (22),  is  only  affected  indi¬ 
rectly  by  representating  x/p  with  | rdrag I Udl  thus,  this  should  not 
act  as  a  missing  sink  of  energy  overall.  The  energy  budget  will  be 
discussed  in  more  detail  in  Section  4.5. 

4.3.  Inline  estimates  of  the  input/output  energy  budget  terms 

In  this  section,  we  evaluate  two  different  inline  estimates  of  the 
input  and  output  terms  in  the  energy  budget,  (22),  one  from  a  sim¬ 
ulation  without  wave  drag  and  one  from  a  simulation  with  wave 
drag.  For  the  simulation  with  wave  drag,  the  C05  scheme  is  used 
in  conjunction  with  the  inline  stratification  and  velocities  from  HY¬ 
COM.  The  four  different  globally  integrated  energy  dissipation 
rates  and  the  wind  power  put  into  the  geostrophic  and  total  flows 
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Fig.  4.  Offline  estimates  of  the  momentum  terms  associated  with  wave  drag,  x/p  (black)  from  (11 )  and  rdm  iid  (red)  from  (14),  in  units  [m2  s  2],  Also  shown  are  the  (panels  a, 
c,  and  e)  differences  between  the  wave  drag  vector  angles  (direction  of  x/p  relative  to  east  minus  direction  of  | rdraK  | ud  relative  to  east)  in  units  of  [°]  and  the  (panels  b,  d,  and  f) 
differences  between  the  wave  drag  vector  magnitudes  (magnitude  of  x/p  minus  magnitude  of  \rdrag\Ud)  in  units  [m2  s  2]  in  the  Southern  Ocean  (a,  b),  equatorial  Pacific  Ocean 
(c,  d),  and  North  Atlantic  Ocean  (e,  f).  Computations  shown  here  were  performed  with  an  example  five-day  average  from  the  last  year  of  the  1  /12°  HYCOM  simulation  without 
wave  drag  and  the  vectors  are  shown  separated  by  ten  degree  spacings. 


Table  2 

The  estimated  globally  integrated  energy  dissipation  rate  [TW  =  1012W]  due  to 
quadratic  bottom  drag  (BD),  parameterized  internal  lee  wave  drag  (WD),  vertical 
viscosity  (W),  and  horizontal  viscosity  (HV)  averaged  over  the  final  year  of  the  spin- 
up  (Year)  of  HYCOM  with  wave  drag  and  the  final  year  of  the  spin-up  (Year)  of 
HYCOM  without  wave  drag.  Also  listed  are  the  wind  power  inputs  to  the  total  flow 
(WindT)  and  geostrophic  flow  (WindC),  the  latter  of  which  excludes  the  regions 
within  five  degrees  of  the  equator. 


Wave  drag? 

Year 

WindT 

WindG 

BD 

WD 

w 

HV 

Yes 

20 

0.868 

0.626 

0.140 

0.402 

0.275 

0.257 

No 

14 

0.867 

0.643 

0.307 

N/A 

0.286 

0.294 

are  listed  in  Table  2.  Maps  of  the  wind  power  put  into  the  geo¬ 
strophic  and  total  flows  are  given  in  Fig.  5a  and  c,  respectively.  Dif¬ 
ferences  between  the  wind  power  input  to  the  total  and 
geostrophic  terms  in  the  simulations  with  versus  without  wave 
drag  are  shown  in  Fig.  5b  and  d,  respectively. 

The  energy  inputs  (Fig.  5a  and  c)  take  on  some  of  their  largest 
values  in  the  western  boundary  current  and  ACC  regions.  The  wind 
power  put  into  the  geostrophic  flow  is  about  0.63  TW,  0.24  TW  less 
than  the  wind  power  input  to  the  total  flow  (Table  2).  Gnanadesi- 
kan  et  al.  (2005)  find  a  similar  difference  between  their 
model-based  estimates  for  the  wind  power  input  to  the  total  and 
geostrophic  flows  -  their  Table  1 .  The  total  surface  velocities  and 
sea  surface  height  distribution  are  locally  altered  with  the  addition 


of  wave  drag,  as  shown  in  Fig.  5b  and  d,  respectively.  Flowever,  the 
global  integrals  of  the  wind  input  into  the  total  and  geostrophic 
flows  are  not  substantially  altered.  It  is  possible  that  the  interan¬ 
nual  variability  in  the  surface  velocities  and  sea  surface  height  dis¬ 
tribution  is  responsible  for  the  differences  we  find  in  Fig.  5b  and  d, 
but  similar  differences  are  found  using  output  from  a  different  year 
of  the  simulation  with  wave  drag  (not  shown). 

Because  we  have  not  computed  the  conversion  from  kinetic  to 
potential  energy9  in  either  simulation,  and  we  have  not  computed 
the  material  derivatives  of  kinetic  and  potential  energy  in  the  simu¬ 
lation  without  wave  drag,  we  cannot  attribute  the  changing  imbal¬ 
ance  between  inputs  and  dissipators  in  (16)  upon  insertion  of 
wave  drag  into  the  model  (Table  2)  to  a  particular  term.  However, 
as  we  noted  in  Section  2,  the  globally  integrated  kinetic  energy  de¬ 
creases  and  the  globally  integrated  potential  energy  increases  upon 
insertion  of  wave  drag.  Also,  the  material  derivative  of  the  mechan¬ 
ical  energy  should  asymptote  to  a  relatively  small  value  as  the  model 
spins  up,  which,  as  we  will  see  in  Section  4.5,  is  the  case  in  the 
simulation  with  wave  drag.  Thus,  it  is  likely  that  the  conversion  from 
kinetic  energy  to  potential  energy  term  in  (16)  compensates  for  the 
increase  in  net  dissipation  upon  insertion  of  wave  drag,  while  the 
wind  power  input  stays  relatively  constant. 


9  We  did  not  compute  the  conversion  from  kinetic  to  potential  energy  due  to 
difficulties  with  computing  the  vertical  velocities,  which  are  not  required  as 
prognostic  variables  in  hybrid  coordinate  models,  below  the  surface  in  HYCOM. 
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(a)  Wind  input  (total  flow)  [W  nr] 
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(c)  Wind  input  (geostrophic  flow)  [W  m 
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(d)  Wind  input  (geostrophic  flow)  percent  difference 


Fig.  5.  The  wind  input:  (a)  to  the  total  flow  and  (c)  to  the  geostrophic  flow,  where  each  power  input  term  has  units  [W  m  Shown  is  an  average  of  inline  estimates  over  the 
final  year  of  the  spin-up  phase  without  wave  drag.  The  percent  differences  between  the  simulations  (ratio  of  the  difference,  without  minus  with  wave  drag,  to  their  sum)  in 
the  wind  input:  (b)  to  the  total  flow  and  (d)  to  the  geostrophic  flow  are  also  shown.  An  average  of  inline  estimates  over  the  final  year  of  the  spin-up  phase  of  each  simulation 
has  been  used.  Listed  in  Table  2  are  the  globally  integrated  power  contributions  of  each  term  shown  in  panels  a  and  c. 


Some  spatial  features  of  the  vertical  and  horizontal  viscosity 
terms  are  shown  in  Fig.  6  and  discussed  below  because  their 
large-scale  patterns  can  be  physically  relevant.  Supporting  the 
hypotheses  of  Thomas  and  Taylor  (2010)  and  Zhai  et  al.  (2012)  that 
a  combination  of  mixed  layer  processes  (e.g.,  shear  generated  in 
the  mixed  layer)  could  leave  less  energy  available  to  deep  ocean 
mixing  than  previously  thought,  we  find  a  relatively  large  dissipa¬ 
tion  rate  due  to  vertical  viscosity  in  the  mixed  layer  alone  (Fig.  6b). 
More  than  three  quarters  (79.7%)  of  the  energy  dissipated  by  ver¬ 
tical  viscosity  is  dissipated  in  the  mixed  layer.  The  zonally-aver- 
aged  viscosity  dissipation  terms  (Fig.  6b  and  d)  in  each  isopycnal 
layer  have  been  multiplied  by  their  respective  layer  thicknesses 
to  make  their  units  consistent  with  the  depth-integrated  viscosity 
dissipation  terms  (Fig.  6a  and  c).  The  vertical  viscosity  term  is  larg¬ 
est  in  equatorial  regions  (Fig.  6a),  in  the  intensified  jets  (Fig.  6a), 
and  in  the  mixed  layer  (Fig.  6b).  The  horizontal  viscosity  term  is 
largest  along  the  fringes  of  the  intensified  jets  (Fig.  6c),  in  the  deep 
convection  and  overflow  regions  (Fig.  6c),  and  along  bathymetric 
boundaries  (Fig.  6d).  The  vertical  viscosity  term  is  always  largest 
in  the  hemisphere  that  is  experiencing  wintertime  conditions 
(which  can  be  inferred  from  Fig.  7)  and  is  the  only  dissipative  term 
that  seasonally  varies  by  an  order  of  0.1  TW,  the  same  order  of 
magnitude  of  the  seasonal  variability  in  the  wind  input  (Fig.  7).  Evi¬ 
dence  of  deep  convection  (and  the  formation  of  overflow  waters  in 


the  northern  hemisphere)  can  clearly  be  seen  in  the  plumes  of 
large  zonally  averaged  vertical  viscosity  dissipation  rates  that  go 
to  about  2000  m  depth  around  50°-55°N  and  about  3000  m  depth 
around  55°-60°S  (Fig.  6b).  Large  zonally  averaged  horizontal  vis¬ 
cosity  dissipation  rates  outline  the  regions  of  deep  convection 
and  formation  of  overflow  waters  (Fig.  6d). 

We  find  that  point  estimates  of  both  of  the  viscosity  terms  can 
be  unphysical  (i.e.,  that  there  are  single  grid  points  that  suggest 
dissipative  processes  are  adding  energy  to  the  system).  Shown  in 
Fig.  6  are  the  negative  of  each  of  the  viscosity  terms,  so  a  positive 
viscosity  dissipation,  for  example,  is  shown  as  a  negative  value  in 
Fig.  6.  At  a  single  grid  point,  there  is  no  physical  meaning  to  a  hor¬ 
izontal  viscosity,  let  alone  energy  dissipation  rate  due  to  horizontal 
viscosity  (Fox-Kemper  and  Menemenlis,  2008).  As  with  the  hori¬ 
zontal  viscosity  term,  the  vertical  viscosity  term  can  be  either  po¬ 
sitive  or  negative  locally,  but  is  negative-definite  when  globally 
integrated.  The  vertical  viscosity  term  is  positive  following  the 
jet  of  the  ACC  near  the  surface  (Fig.  6a),  the  regions  and  depths 
in  the  northern  hemisphere  to  which  water  subducts,  and  in  the 
locations  of  Antarctic  Bottom  Water  (Fig.  6b).  The  horizontal 
viscosity  term  has  some  positive  values  along  the  coasts  and  in 
the  middle  of  the  open  ocean  (Fig.  6c).  Also,  it  can  be  argued  that, 
at  coarser  horizontal  grid  resolutions,  the  larger  mid-depth  hori¬ 
zontal  viscosity  required  leaves  less  energy  to  be  dissipated  in 
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Fig.  6.  The  (a,  b)  vertical  viscosity  and  (c,  d)  "horizontal"  viscosity  terms  in  our  energy  budget,  where  each  viscosity  term  has  units  [W  m  2]:  (a  and  c)  show  depth-integrated 
viscosity  terms  and  (b  and  d)  show  zonally-averaged  viscosity  terms  as  a  function  of  depth.  Shown  is  an  average  of  inline  estimates  over  the  final  year  of  the  spin-up  phase 
without  wave  drag.  Listed  in  Table  2  are  the  globally  integrated  power  contributions  of  each  term.  Shown  with  magenta  contours  (panel  b)  are  the  deepest  mixed  layer 
depths,  as  determined  by  KPP,  from  the  final  year  of  the  1/12°  HYCOM  simulation  without  wave  drag. 


the  bottom,  thus  implying  weaker  abyssal  velocities.  The  same 
argument  can  be  applied  to  coarser  vertical  grid  resolutions,  where 
the  larger  vertical  viscosity  required  near  the  surface  leaves  less 
energy  to  be  dissipated  in  the  abyss,  thus  implying  weaker  abyssal 
velocities.  This  suggests  that  in  simulations  with  higher  grid  reso¬ 
lutions,  our  estimates  of  bottom  drag  and  wave  drag  dissipation 
rates  should  increase,  but  only  up  to  a  point.  This  will  be  investi¬ 
gated  in  a  future  study. 

Differences  between  the  vertical  and  horizontal  viscosity  terms 
with  and  without  wave  drag  are  shown  in  Fig.  8.  In  the  simulation 
with  wave  drag,  the  vertical  viscosity  term  is  smaller  in  regions 
where  deep  convection  occurs  but  larger  in  the  mixed  layer  when 
compared  to  the  simulation  without  wave  drag  (Fig.  8a  and  b).  The 
horizontal  viscosity  term  is  larger  in  regions  where  deep  convec¬ 
tion  occurs  but  smaller  near  the  seafloor  in  the  simulation  with 
wave  drag  than  in  the  simulation  without  wave  drag  (Fig.  8c  and 
d).  However,  in  both  simulations,  the  global  integral  of  the 
horizontal  viscosity  term  is  comparable  in  magnitude  to  that  of 
the  vertical  viscosity  term  (Table  2).  While  Fig.  8  suggests  that 
the  partition  of  energy  is  altered  locally  when  wave  drag  is  in¬ 
serted,  Table  2  suggests  that,  with  the  exception  of  the  bottom  drag 
term,  the  partition  of  energy  dissipated  by  viscosity  is  only  margin¬ 
ally  impacted  globally  when  wave  drag  is  inserted. 

The  differences  between  our  vertical  viscosity  dissipation  rates 
with  and  without  wave  drag  (Fig.  8a  and  b)  suggest  that  adding 


Power  input/output  [TW=1012W] 


Fig.  7.  The  seasonal  cycle  of  the  estimated  energy  dissipation  rate  [TW  =  1012  W] 
due  to  quadratic  bottom  drag,  parameterized  internal  lee  wave  drag,  vertical 
viscosity,  and  horizontal  viscosity,  from  the  last  year  of  a  climatologically  forced 
spin-up  of  HYCOM  with  wave  drag,  as  a  function  of  climatological  month.  Also 

shown  is  the  wind  power  input  to  the  total  flow.  Here,  T’  refers  to  January . and 

'12'  refers  to  December. 
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(a)  Depth-integrated  vertical  viscosity  difference  (without-with  wave  drag) 


(c)  Depth-integrated  horizontal  viscosity  difference  (without-with  wave  drag) 
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Fig.  8.  The  average  difference  in  the  energy  dissipation  rates  by  (a,  b)  vertical  viscosity  and  (c,  d)  “horizontal"  viscosity  terms  in  our  energy  budget  from  the  final  year  of  the 
1/12”  HYCOM  simulation  without  wave  drag  minus  from  the  final  year  of  the  simulation  with  wave  drag,  where  each  panel  has  units  [W  nr2].  Here,  (a  and  c)  show  depth- 
integrated  viscosity  terms  and  (b  and  d)  show  zonally-averaged  viscosity  terms  as  a  function  of  depth.  Positive  values  imply  a  reduction  in  energy  dissipation  and  negative 
values  imply  an  increase  in  energy  dissipation  upon  addition  of  wave  drag  to  the  model.  Shown  with  black  contours  (panel  b)  are  the  deepest  mixed  layer  depths,  as 
determined  by  KPP,  from  the  final  year  of  the  1/12”  HYCOM  simulation  without  wave  drag.  Shown  with  green  contours  (panel  b)  are  the  deepest  mixed  layer  depths,  as 
determined  by  KPP,  from  the  final  year  of  the  1  /12”  HYCOM  simulation  with  wave  drag. 


wave  drag  to  the  momentum  equations  enhances  the  vertical  dif- 
fusivity  in  a  number  of  locations  that  previous  studies  have  exam¬ 
ined.10  For  example,  Polzin  et  al.  (1997)  and  St.  Laurent  et  al.  (2012) 
found  enhanced  levels  of  diapycnal  diffusivity  above  rough  topogra¬ 
phy  in  the  Brazil  Basin  and  in  the  Drake  Passage,  respectively.  Also, 
Ledwell  et  al.  (2000)  and  Kunze  et  al.  (2006)  found  enhanced  levels 
of  energy  dissipation  above  rough  topography  around  the  northern 
portion  of  the  Mid-Atlantic  Ridge  and  just  south  of  Iceland. 

Our  bottom  drag  dissipation  rates  are  more  comparable  to  those 
from  previous  studies  when  we  insert  wave  drag.  Maps  of  the  bot¬ 
tom  drag  and  wave  drag  dissipation  rates  per  unit  area  are  given  in 
Fig.  9a  and  c,  respectively.  Maps  of  the  coefficients  associated  with 
bottom  drag  (Ca\ub\)  and  wave  drag  (rdra g)  are  given  in  Fig.  9b  and  d, 
respectively.  While  we  find  bottom  drag  dissipation  rates  in  the 
simulation  without  wave  drag  (0.31  TW;  see  Table  2)  to  be  larger 
than  a  previous  estimate  of  0.2  TW  based  on  an  observationally-in- 
formed  least-squares  regression  fit  (Sen  et  al.,  2008),  our  energy 
dissipation  rates  due  to  bottom  drag  are  smaller  than  this  previous 
estimate  when  wave  drag  is  inserted  (0.14  TW;  see  Table  2).  Our 


10  Each  component  of  {dvz/dz)(d u/dz)  tends  to  be  at  least  an  order  of  magnitude  less 
than  each  component  of  t ^(cfu/di2)  in  (17)  which  suggests  that  we  may  interpret  the 
relative  magnitude  of  the  vertical  viscosity  term  in  our  energy  budget  as  a  relative 
magnitude  of  the  vertical  diffusivity. 


energy  dissipation  rates  due  to  bottom  drag  in  the  simulation 
without  wave  drag  are  larger  than  the  Arbic  et  al.  (2009)  estimates 
based  on  two  non-data-assimilative  high-resolution  models,  but 
our  bottom  drag  estimates  are  comparable  to  the  estimates  of  Ar¬ 
bic  et  al.  (2009)  when  wave  drag  is  inserted.  Our  energy  dissipation 
rates  due  to  bottom  drag  may  be  larger  than  those  of  Aiki  et  al. 
(2011)  because  Aiki  et  al.  (2011)  defined  bottom  friction  relative 
to  the  average  vertical  mixing  (as  opposed  to  the  quadratic  drag 
used  in  our  model).  The  bottom  drag  and  wave  drag  terms  in  our 
energy  budget  are  large  in  the  intensified  jet  regions,  in  the  over¬ 
flow  regions  around  Iceland,  in  the  Labrador  Sea,  on  the  coastal 
margins  of  the  Arctic  Ocean  north  of  Asia,  around  the  Nordic  Seas, 
and  in  the  Bering  Strait  (Fig.  9a  and  c). 

The  percent  differences  between  the  offline  and  inline  estimates 
of  dissipation  rates  associated  with  wave  drag  are  shown  in  Fig.  10. 
About  half  of  the  dissipation  associated  with  wave  drag  occurs  in 
regions  where  there  are  no  abyssal  hills,  which  increases  the  glob¬ 
ally  integrated  dissipation  rate,  but  the  local  dissipation  rates  due 
to  wave  drag  are  generally  smaller  when  calculated  inline  instead 
of  offline  (Fig.  10).  In  regions  where  there  are  abyssal  hills,  the  in¬ 
line  wave  drag  energy  dissipation  rates  (Fig.  9c)  are  spatially  sim¬ 
ilar  to  those  from  offline  estimates  (Fig.  3b)  with  one  noteworthy 
exception.  The  dissipation  rates  are  substantially  reduced  along 
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(a)  Bottom  drag[login(WmL)] 


(b)  Cd  |uj  [log10(m  s1)] 


(c)  Wave  drag  [log1Q(W  m  ^)] 
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Fig.  9.  The  logio  of  the  bottom  and  wave  drag  terms  in  our  energy  budget,  where  each  drag  term  has  units  [W  m  2]:  (a)  quadratic  bottom  boundary  layer  drag  and  (c) 
parameterized  internal  lee  wave  drag.  Shown  is  an  average  of  inline  estimates  over  the  final  year  of  the  spin-up  phase  with  wave  drag.  Listed  in  Table  2  are  the  globally 
integrated  power  contributions  of  each  term  in  panels  a  and  c.  Also  shown  is  the  log10  of  the  bottom  and  wave  drag  coefficients,  where  each  drag  term  has  units  [m  s-1]:  (b) 
Cd\ub\  and  (d)  rdrag. 
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Fig.  10.  The  average  percent  difference  (ratio  of  the  difference  to  the  sum)  between  the  wave  drag  term  in  our  energy  budget  from  the  offline  estimates  that  make  use  of  the 
5-day  average  velocities  and  stratification  in  the  bottom  500  m  from  the  final  year  of  the  1/12°  HYCOM  simulation  without  wave  drag  versus  from  the  inline  estimates 
computed  during  every  2-min  baroclinic  time  step  during  the  final  year  of  the  1  /12°  HYCOM  simulation  with  wave  drag. 


the  equator  when  wave  drag  is  inserted  into  the  model.  A  compar¬ 
ison  of  the  dissipation  rates  from  Nikurashin  and  Ferrari  (2011) 
versus  those  in  Scott  et  al.  (2011)  and  in  our  offline  estimates  sug¬ 
gests  that  allowing  for  an  anisotropic  spectral  form  of  the  abyssal 


hill  power  spectrum  yields  smaller  dissipation  rates  along  the 
equator  than  use  of  an  isotropic  spectral  form  of  the  abyssal  hill 
power  spectrum.  Fig.  10  suggests  that  putting  wave  drag  into  a 
model  reduces  the  dissipation  rates  by  up  to  an  order  of  magnitude 
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Fig.  11.  The  log10  of  the  (a,  b)  buoyancy  frequency,  which  has  units  [s  ’],  and  (c,  d)  kinetic  energy,  which  has  units  [m2s  2],  in  the  bottom  layer  of  each  1/12°  HYCOM 
simulation,  averaged  over  the  last  year  of  the  spin-up:  (a  and  c)  without  wave  drag  and  (b  and  d)  with  wave  drag. 


in  equatorial  regions.  As  far  as  the  authors  and  E.  Kunze  (personal 
communication)  are  aware,  there  is  no  evidence  for  enhanced  dis¬ 
sipation  near  the  seafloor  at  the  equator  relative  to  surrounding 
latitudes.  Thus,  our  inline  estimates  at  the  equator  may  be  an 
improvement  upon  the  offline  estimates.  The  regions  where  the 
dissipation  rates  are  increased  when  wave  drag  is  inserted  into 
the  model  relative  to  the  offline  estimates  are  in  the  ACC,  some  ex¬ 
treme  latitudes  in  the  Southern  Hemisphere,  and  the  North  Pacific 
Ocean.1 1  The  net  effect  is  that  the  inline  wave  drag  dissipation  rate 
(about  0.4  TW)  has  a  relatively  small  global  integral  compared  to  the 
offline  estimates  (about  1.2  TW  from  the  offline  calculation  over  all 
regions  using  the  HYCOM  stratification). 

Previous  studies  have  noted  the  possibility  that  bottom  drag  may 
substitute  for  wave  drag  in  models  without  a  wave  drag  parameter¬ 
ization  (Arbic  and  Flierl,  2004;  Wright  et  al.,  2013).  The  visually 
apparent  correlation  between  bottom  drag  and  wave  drag  in 
Fig.  9a  and  c  suggests  that  bottom  drag  may  be  a  plausible  substi¬ 
tute  for  wave  drag.  A  least  squares  linear  regression  fit  (not  shown), 
regressing  energy  dissipation  rates  due  to  bottom  drag  on  those  due 
to  wave  drag,  suggests  that  multiplying  the  energy  dissipation  rates 


1 1  One  reason  for  this  is  that,  because  the  velocities  are  not  linearly  related  to  the 
wave  drag,  using  a  five-day  average  of  two-minute  varying  velocities  to  calculate  the 
dissipation  rates  can  lead  to  a  net  bias. 


due  to  bottom  drag  by  a  factor  of  about  five  would  be  a  substitute 
for  putting  an  active  wave  drag  scheme  in  an  ocean  model.  How¬ 
ever,  the  globally  integrated  energy  dissipation  rates  due  to  bottom 
drag  are  about  a  factor  of  about  2.5  times  less  than  the  globally  inte¬ 
grated  energy  dissipation  rates  due  to  wave  drag.  This  is  because 
the  least  squares  linear  regression  fit  is  influenced  primarily  by 
the  locations  where  there  are  large  energy  dissipation  rates  due 
to  both  bottom  drag  and  wave  drag.  Where  energy  dissipation  rates 
due  to  wave  drag  are  below  about  1CT4  W  nT2,  the  approximately 
linear  relationship  no  longer  holds  (not  shown).  At  locations  with 
weak  wave  drag  dissipation,  bottom  drag  can  be  strong  (e.g.,  north¬ 
east  of  Russia  in  Fig.  9a  and  b).  There  are  additional  reasons  for  the 
fact  that  bottom  drag  cannot  serve  as  a  perfect  substitute  for  wave 
drag.  Bottom  drag  dissipation  is  proportional  to  |ub|3,  while  wave 
drag  dissipation  is  proportional  to  |ud|2.  Furthermore,  Cd|Si,|  is  com¬ 
pletely  determined  by  the  velocities  in  the  bottom  HBD  meters, 
while  rdrag  is  determined  from  properties  related  to  the  underlying 
topography,  the  stratification,  and  the  velocities  in  the  bottom  HWD 
meters.  Note  also  that  HBD  #  HWD. 

In  contrast  to  the  100%  local  dissipation  due  to  bottom  drag, 
Waterman  et  al.  (2013)  found  a  mismatch  between  energy  gener¬ 
ation  and  near-bottom  dissipation  and  suggested  that  this  is  be¬ 
cause  about  80%  of  the  energy  is  dissipated  non-locally.  They 
used  observed  abyssal  stratification  and  abyssal  velocities  north 
of  the  Kerguelen  Island  region  to  perform  the  same  offline  analysis 
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(a)  Zonally  averaged  N  without  wave  drag  [log1  Q(s  )]  (b)  Zonally  averaged  N  difference  (without-with  wave  drag) 

-2 


_25  1000 
E 

,  r  2000 

Q. 

0  3000 
-3.5 


4000 


(c)  Zonally  averaged  KE  without  wave  drag  [log10(m2  s  2)](d)  Zonally  averaged  KE  difference  (without-with  wave  drag) 
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Fig.  12.  The  (a)  log10  of  the  zonally  averaged  buoyancy  frequency  (N),  which  has  units  [s  1  ],  from  the  1  / 1 2  HYCOM  simulation  without  wave  drag.  The  (b)  percent  difference 
in  the  zonally  averaged  N  in  the  1  /12°  HYCOM  simulation  without  wave  drag  versus  in  the  simulation  with  wave  drag  (ratio  of  their  difference  to  their  sum).  The  (c)  logi0  of 
the  kinetic  energy  (KE),  which  has  units  [m2  s  2|,  from  the  1  /12"  HYCOM  simulation  without  wave  drag.  The  (d)  percent  difference  in  the  zonally  averaged  KE  in  the  1  /12° 
HYCOM  simulation  without  wave  drag  versus  in  the  simulation  with  wave  drag  (ratio  of  their  difference  to  their  sum).  Shown  with  black  contours  (all  panels)  are  the  zonally 
and  temporally  averaged  depth  of  the  top  of  the  bottom  layer  in  the  final  year  of  the  1  /12°  HYCOM  simulation  without  wave  drag.  Shown  with  green  contours  (panels  b  and 
d)  are  the  zonally  and  temporally  averaged  depth  of  the  top  of  the  bottom  layer  in  the  final  year  of  the  1  /12°  HYCOM  simulation  with  wave  drag. 


as  Nikurashin  and  Ferrari  (2011).  However,  the  isotropic  form  of 
the  abyssal  hill  power  spectra  used  by  Nikurashin  and  Ferrari 
(2011)  and  Waterman  et  al.  (2013)  contrast  with  the  anisotropic 
power  spectra  employed  by  Scott  et  al.  (2011)  and  in  the  present 
manuscript.  Also,  we  find  that  the  inline  energy  dissipation  due 
to  wave  drag  is  reduced  by  an  order  of  magnitude  when  compared 
to  the  offline  estimates  in  the  region  studied  by  Waterman  et  al. 
(2013).  This  raises  the  question  of  whether  discrepancies  between 
offline  energy  dissipation  estimates  and  observed  near-bottom  dis¬ 
sipation  are  entirely  due  to  non-local  energy  dissipation  or  to  more 
subtle  issues  with  how  the  energy  dissipation  rates  are  estimated. 
Non-local  dissipation  of  energy  due  to  wave  drag  will  be  the  sub¬ 
ject  of  investigation  in  a  future  study. 

4.4.  Impact  of  wave  drag  on  the  abyssal  currents  and  stratification 

The  inline  wave  drag  estimates  are  smaller  than  the  fully  global 
offline  wave  drag  estimates  (i.e.,  the  offline  estimates  that  include 
both  abyssal  hill  and  non-abyssal  hill  regions),  in  part,  because  of 
the  active  feedback  that  generally  reduces  the  bottom  stratification 


in  the  model  (Fig.  11a  and  b),  particularly  in  the  Arctic  and  in  regions 
where  wave  drag  is  strong.  Inline  wave  drag  also  reduces  the  bottom 
current  speeds  in  the  model,  as  demonstrated  in  Fig.  11c  and  d.  The 
area-averaged  bottom  buoyancy  frequency  and  kinetic  energy  in 
the  simulation  with  wave  drag  are  reduced  by  59.5%  and  58.5%, 
respectively,  from  those  in  the  simulation  without  wave  drag.  The 
bottom  stratification  (Fig.  12a)  and  kinetic  energy  (Fig.  12c)  are  re¬ 
duced,  on  average,  when  wave  drag  is  inserted,  but  they  increase  in 
some  locations.  The  locations  with  larger  bottom  kinetic  energy  in 
the  simulation  with  wave  drag  are  in  the  Arctic.  The  locations  with 
larger  bottom  stratification  in  the  simulation  with  wave  drag 
(Fig.  11a  and  b)  are  on  the  slopes  of  topographic  features  that  stand 
a  few  kilometers  higher  than  their  base  (Fig.  2c),  particularly  in  the 
Weddell  Sea  and  near  the  Hawaiian  Ridge. 

The  zonally  averaged  buoyancy  frequency  and  kinetic  energy 
are  generally  reduced  closer  to  the  seafloor  and  enhanced  above 
the  seafloor.  The  differences  in  the  zonally  averaged  buoyancy  fre¬ 
quency  and  kinetic  energy  when  wave  drag  is  inserted  are  shown 
in  Fig.  12b  and  d,  respectively.  The  zonally  averaged  buoyancy  fre¬ 
quency  in  the  simulation  with  wave  drag  is  enhanced  by  an  aver- 
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(a)  Offline  wave  drag  percent  difference  (WOA  N,  without-with  wave  drag) 
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(b)  Offline  wave  drag  percent  difference  (HYCOM  N,  without-with  wave  drag) 
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Fig.  13.  The  average  percent  difference  (ratio  of  the  difference  to  the  sum)  between  the  wave  drag  term  in  our  energy  budget  from  the  offline  estimates  that  make  use  of  the 
average  velocities  in  the  bottom  500  m  from  the  final  year  of  the  1  /12°  HYCOM  simulation  without  wave  drag  versus  from  the  offline  estimates  that  make  use  of  the  average 
velocities  in  the  bottom  500  m  from  the  final  year  of  the  1/12°  HYCOM  simulation  with  wave  drag.  Shown  are  the  percent  differences  using  (a)  the  WOA  stratification 
(stratification  held  constant)  and  (b)  the  HYCOM  stratification  from  each  respective  run  (stratification  varied). 

Mechanical  energy  budget  terms  [TW=1012W] 


age  of  1 .29%  from  that  in  the  simulation  without  wave  drag.  The 
zonally  averaged  kinetic  energy  in  the  simulation  with  wave  drag 
is  reduced  by  an  average  of  28.0%  from  that  in  the  simulation  with¬ 
out  wave  drag.  The  relatively  coarse  vertical  resolution  of  the  WOA 
product  makes  it  difficult  to  compare  the  bottom  stratification  in 
the  simulations  with  and  without  wave  drag  (Fig.  11a  and  b)  with 
that  from  the  observations  (Fig.  2a).  Comparison  of  the  model  with 
and  without  wave  drag  to  various  observational  datasets  will  be 
the  subject  of  a  future  study. 

While  the  primary  reason  for  the  overall  reduction  in  bottom  ki¬ 
netic  energy  is  that  wave  drag  acts  to  remove  momentum  in  the 
bottom  500  m,  the  reason  for  the  overall  reduction  in  bottom  strat¬ 
ification  is  less  clear.  A  plausible  hypothesis  for  the  reduction  in 
bottom  stratification  could  be  that  inserting  wave  drag  changes 
the  kinetic  energy,  which  in  turn  changes  the  Richardson  number 
in  the  bottom  500  m.  That  ultimately  enhances  the  vertical  diffu- 
sivity  via  KPP  and  decreases  the  abyssal  stratification.  Evidence 
for  this  is  shown  in  Fig.  8b,  which  generally  shows  a  small  increase 
in  energy  dissipation  associated  with  vertical  viscosity  when  wave 
drag  is  inserted.  Additionally,  a  decrease  in  the  along-isopycnal  dif¬ 
fusion  could  result  in  less  production  of  heavier  waters  via  cabbel- 
ing  or  thermobaricity,  thereby  reducing  the  density  contrast 
between  layers  at  the  bottom.  Evidence  for  this  is  shown  in 
Fig.  8d  near  the  seafloor  where  there  is  a  substantial  decrease  in 
energy  dissipation  associated  with  horizontal  viscosity  when  wave 
drag  is  inserted.  Using  models  that  resolve  sub-mesoscale  motions, 


Fig.  14.  The  seasonal  cycle  of  the  estimated  power  [TW  =  1012  W]  associated  with 
the  partial  time  derivative  of  kinetic  energy  (PE|t ame;  "KE  time"  in  legend),  the 
advective  flux  of  kinetic  energy  (PE|I„ d„;  “KE  adv”),  the  partial  time  derivative  of 
potential  energy  (PEpBme;  “PE  time”),  the  advective  flux  of  potential  energy  (PEp0i„; 
“PE  adv”),  pressure  work  (Pressure!  “pressure”),  buoyancy  diffusion  {P diffusive*  "diffu¬ 
sive"),  conversion  from  internal  energy  to  potential  energy  (CE|_>Ep;  “£|-  >  EP"), 
input  (Pmputi  “input"),  and  output  (P„„t put',  "output”)  from  the  last  year  of  a 
climatologically  forced  spin-up  of  HYCOM  with  wave  drag  as  a  function  of 
climatological  month.  Here,  T  refers  to  January . and  ‘12’  refers  to  December. 
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Nikurashin  et  al.  (2013)  and  Abe  and  Nakamura  (2013)  also  found 
that  diapycnal  mixing  is  enhanced  in  regions  where  there  is  rough 
topography  and  where  lee  waves  break.  Here,  we  do  not  explicitly 
resolve  sub-mesoscale  motions,  but  it  is  reassuring  that  the  dia¬ 
pycnal  mixing  is  enhanced,  leading  to  a  reduction  in  abyssal  strat¬ 
ification,  when  wave  drag  is  inserted  into  the  model. 

The  differences  between  the  offline  wave  drag  dissipation  rates 
calculated  from  the  average  velocities  in  the  bottom  500  m  in  the 
simulation  without  wave  drag  and  those  in  the  simulation  with 
wave  drag  are  given  in  Fig.  13.  The  offline  dissipation  rates  are  esti¬ 
mated  using  WOA  stratification  in  Fig.  13a  and  using  the  stratifica¬ 
tion  from  each  of  the  two  model  simulations  in  Fig.  13b.  The 
relatively  similar  appearances  of  Fig.  1 3a  and  b  suggest  that  the 
velocities  have  a  stronger  impact  on  wave  drag  than  the  stratifica¬ 
tion.  The  velocities  in  the  bottom  500  m  are  altered  in  the  simula¬ 
tion  with  wave  drag  such  that  the  globally  integrated  dissipation 
rates  associated  with  wave  drag  are  cut  by  more  than  half.  Some 
noteworthy  locations  where  the  dissipation  rates  are  largely  re¬ 
duced  include  the  Arctic  Ocean,  the  Indian  Sector  of  the  Southern 
Ocean,  and  in  equatorial  regions.  However,  the  stratification  and 
velocities  in  the  bottom  500  m  are  altered  in  the  simulation  with 


wave  drag  such  that  the  globally  integrated  wave  drag  dissipation 
rates  are  only  cut  in  half.  Regions  where  the  stratification  has  a 
relatively  strong  impact  include  the  Arctic  Ocean,  the  Subpolar 
North  Atlantic,  and  in  the  wake  of  Drake  Passage,  as  evidenced 
by  the  differences  between  Fig.  13a  and  b. 

4.5.  Total  mechanical  energy  budget  analysis 

The  seasonality  of  the  terms  in  our  total  mechanical  energy 
budget,  (22),  are  shown  in  Fig.  14.  The  largest  terms  are  the  source 
and  sink  terms  from  the  momentum  equations,  which  are  highly 
correlated  and  balance  each  other  out  to  within  23.7%  (2.31%)  in 
the  simulation  with  (without)  wave  drag.  The  next  most  dominant 
terms  are  the  advective  flux  of  potential  energy  and  partial  time 
derivative  of  potential  energy,  followed  by  the  conversion  from 
internal  energy  to  potential  energy,  the  latter  of  which  is  highly 
correlated  with  the  source  and  sink  terms  from  the  momentum 
equations.  The  buoyancy  diffusion  term  is  of  the  same  order  of 
magnitude  as  the  conversion  from  internal  energy  to  potential 
energy  term.  The  pressure  work  term  is  relatively  small,  which  is 
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Fig.  is.  The  estimated  power  per  unit  area  [Wm2]  associated  with  (a)  pressure  work  ( Ppressure',  “Pressure  work”  in  title),  (b)  buoyancy  diffusion  ( Pdijfusive ;  “Buoyancy 
diffusion”),  (c)  conversion  from  internal  energy  to  potential  energy  (C£j_>£p;  "IE  to  PE  conversion”),  (d)  advective  flux  of  potential  energy  ( P£pa£/V ;  “PE  conversion"),  (e)  the 
advective  flux  of  kinetic  energy  (PEKadv>  “KE  conversion”),  and  (f)  the  partial  time  derivative  of  the  potential  energy  (P£ptjme;  “PE  time”)  from  the  last  year  of  a  climatologically 
forced  spin-up  of  HYCOM  with  wave  drag. 
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Table  3 

The  estimated  globally  integrated  power  [TW  =  1012  W]  associated  with  the  partial  time  derivative  of  the  kinetic  energy  ( P£|[0me ).  the  partial  time  derivative  of  the  potential  energy 
(P£pti„),  the  advective  flux  of  kinetic  energy  (P£il0, i„),  advective  flux  of  potential  energy  (PEpadv),  pressure  work  (Ppresstlre),  buoyancy  diffusion  conversion  from  internal 

energy  to  potential  energy  (C£|_>£,),  input  f  Pinput ),  and  output  ( Poutput )  averaged  over  the  final  year  of  the  spin-up  of  HYCOM  with  wave  drag. 

PEptime  PEptime  PEpadv  PEpadv  Ppressure  Pdiffiision  f-Ej->Ep  Pinput  Poutput 

0.0126  0.146  0.00284  -0.174  2.81  x  1(T6  0.0309  0.0865  0.868  1.06 


expected  because  the  ocean’s  mean  sea  surface  height  would 
otherwise  be  significantly  varying  along  with  the  Earth’s  geoid. 
The  magnitude  of  the  sum  of  the  partial  time  derivative  of  kinetic 
energy  and  potential  energy  terms  suggests  that  there  is  non-stea- 
dy-state  behavior  in  the  model.  The  degree  to  which  the  model  is 
not  in  steady-state  can  further  be  diagnosed  by  looking  at  Fig.  1, 
which  suggests  that  the  time  derivative  of  the  kinetic  energy  is  still 
changing  but  by  very  little. 

With  the  exception  of  the  partial  time  derivative  of  the  kinetic 
energy  term,  the  spatial  distributions  of  each  of  the  non-input/out- 
put  terms  in  (22)  are  given  in  Fig.  15.  The  advective  flux  of  poten¬ 
tial  energy  and  kinetic  energy  terms  (Fig.  15d  and  e),  pressure  work 
term  (Fig.  15a),  and  partial  time  derivative  of  potential  energy  term 
(Fig.  15f)  each  take  on  their  largest  values  near  the  equator.  The 
conversion  from  internal  energy  to  potential  energy  term 
(Fig.  15c)  takes  on  its  largest  values  near  the  intensified  jets.  The 
buoyancy  diffusion  term  (Fig.  15b)  takes  on  its  largest  values  just 
equatorward  of  where  deep  convection  occurs.  The  sum  of  the  con¬ 
version  from  internal  energy  to  potential  energy  and  buoyancy  dif¬ 
fusion  terms  is  very  similar  in  spatial  distribution  to  the  buoyancy 
flux  (not  shown)  because  each  of  these  terms  is  closely  related  to 
the  stratification  via  a  vertical  density  gradient.  The  sum  of  the 
advective  flux  of  potential  energy  term  (Fig.  15d)  and  and  the  par¬ 
tial  time  derivative  of  potential  energy  term  (Fig.  15f)  roughly  bal¬ 
ances  the  sum  of  the  rest  of  the  terms. 

Our  total  mechanical  energy  budget,  (22),  is  approximately  bal¬ 
anced  (Table  3)  with  a  <  10%  imbalance  relative  to  the  sum  of  the 


dissipators.  There  are  several  potential  reasons  for  the  imbalance 
between  the  sum  of  the  partial  time  derivative  and  advective  flux 
terms  (-0.0126  TW)  and  the  sum  of  the  pressure  work,  buoyancy 
diffusion,  internal  to  potential  energy  conversion,  input,  and  out¬ 
put  terms  (-0.0746  TW). 

(1)  While  the  Robert- Asselin  time  filter  suppresses  a 
spurious  time-splitting  computational  mode  associated 
with  the  use  of  the  leap-frog  time-stepping  scheme 
(Griffies  et  al.,  2000),  it  is  an  unmeasured  source  of  dissi¬ 
pation.  An  associated  imbalance  in  surface  forcing  gradu¬ 
ally  leads  to  a  small  violation  of  tracer  (e.g.,  heat) 
conservation  over  time  (Leclair  and  Madec,  2009).  These 
factors  could  account  for  some  of  the  disagreement  in 
our  total  mechanical  energy  budget.  The  dissipation  due 
to  the  Robert-Asselin  filter  in  a  one-layer  M2  tide  run 
of  HYCOM  was  quantified  by  doubling  the  filter  weights. 
It  was  found  that  the  dissipation  is  less  than  0.05  TW,  out 
of  a  total  M2  dissipation  rate  of  about  2.4  TW  (not 
shown). 

(2)  We  used  centered  differences  for  all  of  the  spatial  deriva¬ 
tives  and  our  model  frequently  implements  an  interpola¬ 
tion  scheme  between  density  and  z-coordinates.  These 
numerical  errors  are  propagated  and  may  be  significant 
in  a  non-uniformly  spaced  coordinate  system  such  as  that 
of  a  time-varying  isopycnic  coordinate  model  with  a  tri¬ 
pole  grid. 


(a)  Info,  tensor  (1,1)  GAM  residuals 


(b)  Info,  tensor  (1,2)  GAM  residuals 


(c)  Info,  tensor  (2,1)  GAM  residuals  (d)lnfo.  tensor  (2,2)  GAM  residuals 


-150  -100  -50  0  50  100  150  -150  -100  -50  0  50  100  150 

Longitude  Longitude 


Fig.  16.  The  residuals  of  the  GAM  fits  plotted  against  longitude  for  each  of  the  four  G05  wave  drag  scheme  parameters:  (a-d)  the  four  information  tensor  elements  of  T 
[kgm-2  s  ’]. 
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(a)  (1,1)  component  of  info,  tensor  ( b)  10x(1,2)  component  of  info,  tensor 
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Fig.  17.  The  estimates  of  (a-d)  the  four  information  tensor  elements  of  T  [kg  m  2  s  ']  in  abyssal  hill  regions  from  (10)  and  in  non-abyssal  hill  regions  from  using  a  GAM,  as 
described  in  Appendix  B. 


(3)  HYCOM  consistently  makes  the  approximation  that  the  non- 
steric  sea  surface  height  (i.e.,  bottom  pressure  anomaly)  is  a 
small  fraction  of  the  total  depth.  We  are  using  the  same 
approximation  in  the  total  mechanical  energy  budget,  which 
might  be  more  accurate  if  we  instead  used  the  total  depth 
including  the  non-steric  sea  surface  height. 

(4)  We  may  need  to  explicitly  account  for  the  compressibility  of 
water  from  thermobaricity.  However,  this  term  is  a  global 
integral  over  p(dp/dz)/ p,  which  should  be  small  because 
this  is  of  a  comparable  order  of  magnitude  as  the  buoyancy 
diffusion  term  we  have  computed  here. 

(5)  The  power  associated  with  the  along-isopycnal  diffusion  of 
buoyancy  may  be  non-negligible.  Using  typical  values  of 
each  of  the  variables  in  gV  ■  ( KzVp ),  we  estimate  that  the 
diapycnal  contribution  to  the  diffusion  of  buoyancy  term  is 
of  the  same  order  of  magnitude  as  the  along-isopycnal  con¬ 
tribution  to  the  diffusion  of  buoyancy  term.  This  could 
potentially  account  for  our  <10%  energy  budget  imbalance 
and  is  of  the  correct  sign. 

5.  Conclusions 

The  net  dissipative  and  input  terms  in  a  total  mechanical  energy 
budget  have  been  presented  here  for  an  eddying  global  ocean  mod¬ 
el  with  and  without  an  inline  parameterization  of  topographic 
internal  lee  wave  drag.  Topographic  internal  lee  wave  drag  and 
quadratic  bottom  boundary  layer  drag  each  contribute  a  significant 
amount  of  the  total  energy  dissipation  in  the  model.  Horizontal 
and  vertical  eddy  viscosity  are  both  responsible  for  energy  dissipa¬ 
tion  rates  greater  than  those  associated  with  bottom  drag,  but  less 


than  those  associated  with  wave  drag.  Most  of  the  energy  dissipa¬ 
tion  associated  with  horizontal  viscosity  occurs  at  the  margins  of 
the  intensified  currents  of  the  ocean.  The  energy  dissipation 
associated  with  vertical  viscosity  takes  place  primarily  within  the 
mixed  layer.  The  energy  dissipation  rate  by  bottom  drag  is  reduced 
to  a  greater  extent  upon  addition  of  wave  drag  to  the  model  than  is 
the  energy  dissipation  rate  by  either  of  the  viscosity  terms.  How¬ 
ever,  utilization  of  a  boosted  bottom  drag  cannot  serve  as  a  substi¬ 
tute  for  a  wave  drag  parameterization.  Energy  dissipation  by  wave 
drag  more  than  compensates  for  the  decrease  in  energy  dissipation 
rates  via  bottom  drag,  vertical  viscosity,  and  horizontal  viscosity. 

Our  inline  energy  budget  analysis  is  preceded  by  our  own  offline 
analysis,  which  is  much  faster  to  perform.  In  contrast  to  previous 
offline  estimates  of  globally  integrated  dissipation  by  lee  wave  drag 
(Nikurashin  and  Ferrari,  2011;  Scott  et  al.,  2011),  our  own  offline 
estimates  have  used  two  different  wave  drag  schemes  for  compar¬ 
ison.  Offline  estimates  show  that  the  C05  wave  drag  scheme  yields 
comparable  spatial  dissipation  maps  and  comparable  globally  inte¬ 
grated  dissipation  rates  to  the  B75  scheme  in  regions  where  there 
are  abyssal  hills.  Our  offline  globally  integrated  dissipation  rates 
are  closer  to  those  of  Scott  et  al.  (201 1 )  than  to  those  of  Nikurashin 
and  Ferrari  (2011).  When  extended  to  regions  without  abyssal  hills, 
the  offline  globally  integrated  dissipation  rates  are  approximately 
tripled.  Another  important  contrast  between  our  study  and  previ¬ 
ous  studies  is  that  we  estimate  the  dissipation  rates  inline  as  the 
model  runs  as  well  as  offline.  The  extra  roughness  we  choose  to  in¬ 
sert  in  regions  where  there  are  no  abyssal  hills  doubles  the  inline 
globally  integrated  dissipation  rates.  However,  the  inline  feedback 
that  reduces  the  abyssal  velocities  and  stratification  lowers  the 
globally  integrated  dissipation  rates.  In  the  end,  the  two  factors 
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roughly  cancel  each  other  so  that  our  inline  estimates  of  globally 
integrated  dissipation  rates  by  wave  drag  are  comparable  to  global 
offline  estimates  made  in  previous  studies  (Nikurashin  and  Ferrari, 
2011;  Scott  et  al.,  201 1 ).  The  equator  is  the  region  with  the  largest 
absolute  reduction  in  energy  dissipation  rates  due  to  wave  drag  tak¬ 
ing  place  in  inline  versus  offline  estimates.  The  regions  having  the 
largest  percent  reductions  in  wave  drag  energy  dissipation  rates 
in  inline  versus  offline  estimates  are  generally  regions  where  dissi¬ 
pation  is  relatively  small  (mostly  non-abyssal  hill  regions). 

We  have  evaluated  an  approximately  balanced  total  mechanical 
energy  budget.  There  are  a  number  of  confounding  factors  that 
could  be  responsible  for  the  imbalance  of  our  total  mechanical  en¬ 
ergy  budget.  It  is  possible  that  the  compressibility  effects  are  sig¬ 
nificant  and/or  that  the  along-isopycnal  contributions  to  the 
diffusion  of  buoyancy  term  are  important.  Errors  in  the  numerical 
methods,  and  approximations  used  in  the  model,  may  contribute 
to  the  total  mechanical  energy  budget's  imbalance  as  well. 

We  recommend  several  improvements  to  our  implementation 
of  wave  drag.  First,  the  range  of  relevant  wavenumbers  for  the 
internal  waves  to  not  be  evanescent  can  be  different  for  different 
regions  of  the  ocean,  in  contrast  to  the  fixed  range  we  specified 
by  (6).  Second,  a  more  physically  motivated  estimation  of  the 
roughness  parameters  in  non-abyssal  hill  regions,  taken  from  a  for¬ 
mulation  that  characterizes  the  statistical  features  with  a  power 
spectrum,  would  be  an  improvement  over  the  machine  learning 
algorithm  we  used  here  (Appendix  B).  Third,  the  form  of  the  wave 
drag  parameterization  in  the  momentum  equations  could  be  im¬ 
proved  upon.  The  full  tensor  wave  drag  in  (11)  is  preferable  to 
the  reduced  scalar  form  in  (14)  used  here.  Fourth,  the  assumption 
of  100%  local  dissipation  in  the  bottom  500  m  can  be  relaxed  by 
using  the  depth-dependent  momentum  deposition  procedure 
developed  in  G05.  Fifth,  it  should  be  investigated  how  to  imple¬ 
ment  wave  drag  in  a  subset  of  regions  (e.g.,  only  in  regions  where 
the  seafloor  is  deeper  than  500  m,  only  in  abyssal  hill  regions,  only 
where  f  <  N,  or  only  where  topographic  slopes  are  not  supercriti¬ 
cal)  while  maintaining  numerical  stability.  Sixth,  the  diapycnal  dif- 
fusivity  could  be  adjusted  (e.g.,  Polzin,  2009)  according  to  the 
energy  dissipation  arising  from  the  wave  drag  scheme  in  addition 
to  the  indirect  effects  on  the  diapycnal  diffusivity  arising  from 
changes  in  the  vertical  shear  which  then  activate  changes  through 
KPP.  Lastly,  the  implementation  of  the  C05  scheme  would  benefit 
from  the  development  of  a  procedure  to  constrain  all  of  its  free 
parameters  that  are  directly  related  to  the  magnitude  of  the  wave 
drag.  Alternatively,  a  different  wave  drag  scheme  could  be  used, 
but  different  problems,  in  addition  to  some  of  the  aforementioned 
ones,  would  then  arise. 

Here,  we  have  focused  on  the  impact  of  parameterized  topo¬ 
graphic  internal  lee  wave  drag  on  model  energetics.  We  still  need 
to  examine  whether  the  inclusion  of  this  parameterization  im¬ 
proves  the  performance  of  HYCOM  relative  to  various  in  situ  and 
satellite  observations.  Furthermore,  we  need  to  examine  how  the 
role  of  wave  drag  in  the  energy  budget  changes  as  the  horizontal 
resolution  of  the  model  is  increased.  These  currently  unresolved  is¬ 
sues  will  be  pursued  in  a  follow-up  study. 
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Appendix  A.  Wave  Drag  Parameter  Estimation 

In  order  to  implement  the  wave  drag  scheme  of  Garner  (2005; 
G05  hereafter),  we  need  to  estimate  two  parameters  that  are  as¬ 
sumed  to  be  spatially  homogeneous,  y  and  e,  and  one  parameter 
that  we  allow  to  vary  in  space,  /L  The  first  of  these  relates  the  max¬ 
imum  height  of  topographic  features  to  their  horizontal  scale  and 
the  second  determines  the  number  density  of  features  as  a  function 
of  their  maximum  height.  Using  the  same  method  (linear  regres¬ 
sion  using  functions  of  the  velocity  potential  for  linear  mountain 
waves)  described  in  the  Appendix  of  G05  but  applied  to  the  syn¬ 
thetic  topography  described  in  Section  3.2  generated  on  a  2-min 
grid,  we  find  that  y  =  0.36,  which  is  close  to  the  value  of  0.4  in 
G05.  The  other  uniform  parameter,  e,  is  estimated  as  the  least- 
squares  slope  of  log(n(/i))  versus  log(Ali),  where  n(h)  is  the  number 
of  topographic  height  maxima  in  the  interval  Ah  sampled  from 
10°  x  10°  boxes.  We  find  that  e  =  0.02.  The  examples  in  G05  as¬ 
sumed  e  =  0. 

The  third  parameter,  /),  establishes  the  shape  of  the  vertical 
cross-section  of  the  hills  for  the  dimensional  analysis.  We  keep 
the  mean  value  of  p  =  0.5  used  by  G05  (midway  between  a  trian¬ 
gular  and  flat  summit)  and  focus  on  the  variability  about  the  mean. 
Preliminary  model  runs  suggest  that  the  standard  deviation  should 
be  on  the  order  of  0.25.' :  We  impose  this  variability  by  first  exam¬ 
ining  the  distribution  of  A jSy  =  |hx.+1.  -  hx.  tJ\  +  \hy,J+1  -  hyiJ  t  |  where 
Kj  =  hi+  ij  -  hi  ,j.  hyij  =  hij+1  -  hjj  ,,  and  hu  is  the  gridded  high-res- 
olution  topography  with  i  and  j  denoting  longitude  and  latitude  indi¬ 
ces.  A/1  is  a  measure  of  the  steepness  of  the  local  terrain.  Then  we  set 


jSy  =  0.4  +  0.0006A&J,  (A.l) 


which  has  the  desired  mean  and  standard  deviation. 

There  are  other  tunable  parameters  associated  with  the  wave 
drag  scheme  of  G05.  For  example,  a0  is  a  coefficient  for  the  propa¬ 
gating  component  of  wave  drag  and  ai  is  a  coefficient  for  the 
non-propagating  component  of  wave  drag.  We  follow  G05  and  con¬ 
strain  these  remaining  parameters  to  satisfy  the  relationship, 
(di/ao)  =  (Hcrit/Hcrit.r)  (where  Hmt  and  Hoir.r  are  the  non-dimen- 
sionalized  critical  height  and  its  reference).  G05  assumed  that 
Hcr^r  =  1/9  to  yield  a  maximum  total  drag  of  about  2D*,  where 
D*  is  the  base  flux  of  momentum  with  Han  =  0.7  and  a0  =  1.  We 
follow  suit  here. 


12  Using  the  WOA  bottom  stratification  and  bottom  velocities  from  June  and 
December  of  the  final  year  of  the  HYCOM  simulation  without  wave  drag,  we  found 
that  using  values  for  (I  anywhere  between  0  and  1  yielded  energy  dissipation  rates  to 
be  well  within  a  factor  of  2  of  those  from  Scott  et  al.,  2011. 
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Appendix  B.  Extending  inputs  for  the  G05  scheme  to  regions 
without  abyssal  hills 

The  Goff  and  Arbic  (2010)  and  Goff  (2010)  fields  are  only 
defined  in  regions  where  abyssal  hills  are  defined  according  to  Goff 
and  Jordan  (1988).  However,  the  numerics  of  HYCOM  are  unstable 
if  we  apply  large  wave  drag  at  one  grid  point  and  no  wave  drag  at 
an  adjacent  grid  point.  Also,  there  exists  observational  evidence  (St. 
Laurent  et  al.,  2012)  that  the  roughness  in  some  non-abyssal  hill 
regions  may  be  sufficient  to  generate  lee  waves.  Therefore,  we  pur¬ 
sue  a  method  to  estimate  T  such  that  each  element  of  T  smoothly 
varies  toward  zero  away  from  the  abyssal  hills.  We  use  a  machine¬ 
learning  algorithm  called  a  Generalized  Additive  Model  [GAM; 
using  the  mgcv  package’s  bam  in  R  (Wood,  2006)]  with  the  func¬ 
tional  dependence, 

hj  =  /,,y(M  +/2,u(M)  +/3,u(/U).  (B.l) 

Here,  h  is  the  root  mean  square  topographic  height;  p  is  from  (A.l); 
X  is  longitude;  and /„  fj  for  n  =  1 , 2, 3,  i  =  1 , 2,  and  j  =1,2  are  smooth 
but  nonlinear  functions  determined  by  the  best  fit  subject  to  a  con¬ 
dition  that  the  integral  over  the  square  of  their  second  derivative,  a 
measure  of  “wiggliness"  as  described  in  another  context  by  Tross¬ 
man  et  al.  (2011),  is  not  too  large.  Here,  T,  i ,  T12;  T2  i ,  and  T2,2  are 
elements  of  T.  The  estimates  of  each  T q  are  nearly  insensitive  to  lat¬ 
itudinal  dependencies  because  topographic  roughness  varies  more 
with  longitude  than  latitude.  Eight  GAMs  are  trained  where  there 
are  abyssal  hills  from  Goff  and  Arbic  (2010).  Two  GAMs  are  trained 
for  each  of  the  four  T  if  one  for  regions  east  of  the  Prime  Meridian 
and  one  for  regions  west  of  the  Prime  Meridian.  If  each  Ty  is  esti¬ 
mated  using  only  one  GAM,  then  the  values  of  each  Ty  are  larger 
at  every  location  where  there  are  no  abyssal  hills.  Each  GAM,  with 
residuals  shown  to  be  uniformly  scattered  as  a  function  of  longitude 
in  Fig.  16,  is  used  to  predict  what  the  values  of  the  four  T, j  are 
everywhere  outside  of  the  abyssal  hill  regions.  Fig.  16  also  shows 
that  there  is  no  systematic  bias  in  our  GAM-based  estimates  and 
that  the  standard  error  always  stays  less  than  a  few  tens  of  percent 
of  the  field. 

GAMs  are  designed  to  be  able  to  predict  outside  of  their  training 
data,  which  we  do  here.  We  tend  to  predict  each  Tjj  outside  of  the 
abyssal  hill  regions  using  values  of  at  least  two  of  /I,  h,  and  X  found 
in  regions  where  there  are  abyssal  hills.  For  example,  both  inside 
and  outside  of  the  abyssal  hill  regions  we  can  have  small  values 
for  p  at  a  given  longitude.  However,  the  magnitude  of  h  can  be  dif¬ 
ferent  in  abyssal  hill  regions  from  h  in  non-abyssal  hill  regions 
along  the  same  longitude  near  a  coastal  margin.  The  regions  where 
we  predict  each  Tjj  tend  to  have  small  p  and  widely  varying  h  for 
each  X. 

Lastly,  we  smooth  the  resulting  horizontal  gradients  of  each  of 
the  T,j  along  the  boundaries  of  the  abyssal  hill  regions.  To  do  this, 
we  apply  a  Gaussian  Blackman  filter  (with  f  special,  fwind.2,  and 
filter2  in  Matlab)  using  31  grid  points  everywhere,  zeroing  out 
the  normalized  frequencies  greater  than  0.5  and  less  than  0.1. 
The  original  values  of  each  of  the  Tjj  in  regions  where  there  are 
abyssal  hills  are  reassigned  so  that  each  of  the  Tjj  is  left  unaltered 
by  the  Blackman  filter  in  regions  where  there  are  abyssal  hills. 

Each  element  of  the  information  tensor,  T,  is  shown  in  Fig.  17. 
As  expected,  the  diagonal  components  of  T  (Fig.  17a  and  d)  are  of 
the  same  order  of  magnitude,  and  the  off-diagonal  components 
of  T  (Fig.  17b  and  c)  are  an  order  of  magnitude  less  than  the  diag¬ 
onal  components.  In  general,  the  values  of  T  are  larger  in  the  Atlan¬ 
tic  than  in  the  Pacific.  The  GAMs  find  that  there  are  several 
topographic  features  that  could  be  involved  in  lee  wave  generation 
in  the  North  Atlantic  Ocean  (i.e.,  in  the  areas  surrounding  Green¬ 


land  and  Iceland).  The  values  away  from  the  abyssal  hills  in 
Fig.  1 7  decay  to  zero  within  a  few  degrees  except  around  Green¬ 
land  and  Iceland,  where  there  may  be  enough  roughness  to  gener¬ 
ate  internal  lee  waves. 
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